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Abstract 

We consider the following basic learning task: given independent draws from an unknown 
distribution over a discrete support, output an approximation of the distribution that is as ac¬ 
curate as possible in i\ distance (equivalently, total variation distance, or “statistical distance”). 
Perhaps surprisingly, it is often possible to “de-noise” the empirical distribution of the samples 
to return an approximation of the true distribution that is significantly more accurate than the 
empirical distribution, without relying on any prior assumptions on the distribution. We present 
an instance optimal learning algorithm which optimally performs this de-noising for every dis¬ 
tribution for which such a de-noising is possible. More formally, given n independent draws from 
a distribution p, our algorithm returns a labelled vector whose expected distance from p is equal 
to the minimum possible expected error that could be obtained by any algorithm that knows 
the true unlabeled vector of probabilities of distribution p and simply needs to assign labels, 
up to an additive subconstant term that is independent of p and goes to zero as n gets large. 
This somewhat surprising result has several conceptual implications, including the fact that, for 
any large sample, Bayesian assumptions on the “shape” or bounds on the tail probabilities of a 
distribution over discrete support are not helpful for the task of learning the distribution. 

As a consequence of our techniques, we also show that given a set of n samples from an 
arbitrary distribution, one can accurately estimate the expected number of distinct elements 
that will be observed in a sample of any size up to n log n. This sort of extrapolation is practically 
relevant, particularly to domains such as genomics where it is important to understand how much 
more might be discovered given larger sample sizes, and we are optimistic that our approach is 
practically viable. 


‘This work is supported in part by NSF CAREER Award CCF-1351108. 
^This work is supported in part by a Sloan Research Fellowship. 



1 Introduction 


Given independent draws from an unknown distribution over an unknown discrete support, what 
is the best way to aggregate those samples into an approximation of the true distribution? This 
is, perhaps, the most fundamental learning problem. The most obvious and most widely employed 
approach is to simply output the empirical distribution of the sample. To what extent can one 
improve over this naive approach? To what extent can one “de-noise” the empirical distribution, 
without relying on any assumptions on the structure of the underlying distribution? 

Perhaps surprisingly, there are many settings in which de-noising can be done without a priori 
assumptions on the distribution. We begin by presenting two motivating examples illustrating 
rather different settings in which significant de-noising of the empirical distribution is possible. 

Example 1. Suppose you are given 100,000 independent draws from some unknown distribution, 
and you find that there are roughly 1,000 distinct elements, each of which appears roughly 100 
times. Furthermore, suppose you compute the variance in the number of times the different domain 
elements occur, and it is close to 100. Based on these samples, you can confidently deduce that the 
true distribution is very close to a uniform distribution over 1,000 domain elements, and that the 
true probability of a domain element seen 90 times is roughly the same as that of an element observed 
110 times. The basic reasoning is as follows: if the true distribution were the uniform distribution, 
then the noise from the random sampling would exhibit the observed variance in the number of 
occurrences; if there was any significant variation in the true probabilities of the different domain 
elements, then, combined with the noise added via the random sampling, the observed variance 
would be noticeably larger than 100. The l\ error of the empirical distribution would be roughly 0.1, 
whereas the u de-noised v distribution would have error less than 0.01. 

Example 2. Suppose you are given 1,000 independent draws from an unknown distribution, and all 
1000 samples are unique domain elements. You can safely conclude that the combined probability of 
all the observed domain elements is likely to be much less than 1/100—if this were not the case, one 
would expect at least one of the observed elements to occur twice in the sample. Hence the empirical 
distribution of the samples is likely to have i\ distance nearly 2 from the true distribution, whereas 
this reasoning would suggest that one should ascribe a total probability mass of at most 1/100 to 
the observed domain elements. 

In both of the above examples, the key to the “de-noising” was the realization that the true 
distributions possessed some structure—structure that was both easily deduced from the samples, 
and structure that, once known, could then be leveraged to de-noise the empirical distribution. 
Our main result is an algorithm which de-noises the empirical distribution as much as is possible, 
whenever such denoising is possible. Specifically, our algorithm achieves, up to a subconstant term 
that vanishes as the sample size increases, the best error that any algorithm could achieve—even 
an algorithm that is given the unlabeled vector of true probabilities and simply needs to correctly 
label the probabilities. 

Theorem 1. There is a function err(n) that goes to zero as n gets large, and an algorithm, which 
given n independent draws from any distribution p of discrete support, outputs a labelled vector q, 
such that 

E[\\p — < 7 |1 1 ] < opt(p,n) + err(n), 
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where opt(p, n ) is the minimum expected error that any algorithm could achieve on the following 
learning task: given p, and given n samples drawn independently from a distribution that is identical 
to p up to an arbitrary relabeling of the domain elements, learn the distribution. 

The performance guarantees of the above algorithm can be equivalently stated as follows: let 
S <— p denote that S is a set of n independent draws from distribution p, and let n(p) denote a 

n 

distribution that is identical to p, up to relabeling the domain elements with arbitrary distinct new 
labels given by the mapping it. Our algorithm, which maps a set of samples S to a labelled vector 
q = f(S), satisfies the following: For any distribution p, 

E [||p - g||i] < min max ( E [t r(p) - A(S)} ) + o n (l), 

algs A n \S<—ir(p) I 

where o n (l) —>• 0 as n —>• oo is independent of p and depends only on n. 

One surprising implication of the above result is that, for large samples, prior knowledge of 
the “shape” of the distribution, or knowledge of the rate of decay of the tails of the distribution, 
cannot improve the accuracy of the learning task. For example, typical Bayesian assumptions 
that the frequency of words in natural language satisfy Zipf distributions, or the frequencies of 
different species of bacteria in the human gut satisfy Gamma distributions or various power-law 
distributions, etc, can improve the expected error of the learned distribution by at most a vanishing 
function of the sample size. 

The key intuition behind this optimal de-noising, and the core of our algorithm, is the ability 
to very accurately approximate the unlabeled vector of probabilities of the true distribution, given 
access to independent samples. In some sense, our result can be interpreted as the following 
statement: up to an additive subconstant factor, one can always recover an approximation of the 
unlabeled vector of probabilities more accurately than one can disambiguate and label such a vector. 
That is, if one has enough samples to accurately label the unlabeled vector of probabilities, then 
one also has more than enough samples to accurately learn that unlabeled vector. Of course, this 
statement can only hold up to some additive error term, as the following example illustrates. 

Example 3. Given samples drawn from a distribution supported on two unknown domain elements, 
if one is told that both probabilities are exactly 1/2, then as soon as one observes both domain 
elements, one knows the distribution exactly, and thus the expected error given n samples will be 
0(l/2 n ) as this bounds the probability that one of the two domain elements is not observed in a 
set of n samples. Without the prior knowledge that the two probabilities are 1/2, the best algorithm 
will have expected error ~ 1 /y/n. 

The above example illustrates that prior knowledge of the vector of probabilities can be helpful. 
Our result, however, shows that this phenomena only occurs to a significant extent for very small 
sample sizes; for larger samples, no distribution exists for which prior knowledge of the vector of 
probabilities improves the expected error of a learning algorithm beyond a universal subconstant 
additive term that goes to zero as a function of the sample size. 

1.1 Our Approach 

Our algorithm proceeds via two steps. In the first step, the samples are used to output an ap¬ 
proximation of the vector of true probabilities. We show that, with high probability over the 
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randomness of the n independent draws from the distribution, we accurately recover the portion 
of the vector of true probabilities consisting of values asymptotically larger than l/(nlogn). Note 
that the empirical distribution accurately estimates probabilities down to ~ 1/ n —indeed the vector 
of empirical probabilities are all multiples of 1/n. The characterization of the first phase of our 
algorithm can be interpreted as showing that we recover the vector of probabilities essentially to 
the accuracy that the empirical distribution would have if it were based on nlogn samples, rather 
than only n samples. Of course, this surprisingly accurate reconstruction comes with the caveat 
that we are only recovering the unlabeled vector of probabilities—we do not know which domain 
elements correspond to the various probabilities. 

The second step of our algorithm leverages the accurate approximation of the unlabeled vector 
of probabilities to optimally assign probability values to each of the observed domain elements. For 
some intuition into this step, first suppose you know the exact vector of unlabelled probabilities. 
Consider the following optimization problem: given n independent draws from distribution p , 
and an unlabeled vector v representing the true vector of probabilities of distribution p, for each 
observed domain element x, assign the probability q(x ) that minimizes the expected t\ distance 
| q(x) — p(x)\. As the following example illustrates, this problem is more subtle than it might initially 
seem; intuitive schemes such as assigning the ith largest probability in v to the element with the 
ith largest empirical probability is not optimal. 

Example 4. Consider n independent draws from a distribution in which 90% of the domain el¬ 
ements occur with probability 10 /n, and the remaining 10% occur with probability 11/n. If one 
assigns probability 11 /n to the 10% of the domain elements with largest empirical frequencies, the 
l\ distance will be roughly 0.2, because the vast majority of the elements with the largest empirical 
frequencies will actually have true probability 10/n rather than 11/n. In contrast, if one ignores the 
samples and simply assigns probability 10/n to all the domain elements, the t\ error will be exactly 
0.1. 


This optimization task of assigning probabilities q(x) (as a function of the true probabilities v 
and set of n samples) so as to minimize the expected I\ error is well-defined. Nevertheless, this task 
seems to be computationally intractable. Part of the computational challenge is that the optimal 
probability to assign to a domain element x might be a function of v, the number of occurrences of 
x in the sample, and also the number of occurrences of all the other domain elements. Nevertheless, 
we describe a very natural and computationally efficient scheme which assigns a probability q(x) 
to each x that is a function of only v and the number of occurrences of x, and we show that this 
scheme incurs an expected error within o(l) of the expected error of the optimal scheme (which 
assigns q(x) as a function of v and the entire set of samples). Of course, there is the additional 
complication that, in the context of our two step algorithm, we do not actually have the exact 
vector of true probabilities—only an approximation of such a vector—and hence this second phase 
of our algorithm must be robust to some noise in the recovered probabilities. 

Beyond yielding a near optimal learning algorithm, there are several additional benefits to our 
approach of first accurately reconstructing the unlabeled vector of probabilities. For instance, 
such an unlabeled vector allows us to estimate properties of the underlying distribution including 
estimating the error of our returned vector, and estimating the error in our estimate of each observed 
domain element’s probability. Additionally, as the following proposition quantifies, this unlabeled 
vector of probabilities can be leveraged to produce an accurate estimate of the expected number of 
distinct elements that will be observed in sample sizes up to a logarithmic factor larger: 
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Proposition 1. Given n samples from an arbitrary distribution p, with high probability over the 
randomness of the samples, one can estimate the expected number of unique elements that would be 
seen in a set of k samples drawn from p, to within error k ■ cJ n £ n for some universal constant c. 


This proposition is tight, and it is slightly surprising in that the factor by which one can 
accurately extrapolate increases with the sample size. This ability to accurately predict the expected 
number of new elements observed in larger sample sizes is especially applicable to such settings as 
genomics, where data is relatively costly to gather, and the benefit of data acquisition is largely 
dependent on the number of new phenomena discovered^] 


1.2 Related Work 

Perhaps the first work on correcting the empirical distribution is the work of Turing, and I.J. 
Good [21j (see also [22] ) —which serves as the jumping-off point for nearly all of the subsequent 
work on this problem that we are aware. In the context of their work at Bletchley Park as part of 
the British WWII effort to crack the German enigma machine ciphers, Turing and Good developed 
a simple estimator that corrected the empirical distribution, in some sense to capture the “missing” 
probability mass of the distribution. This estimator and its variants have been employed widely, 
particularly in the contexts of genomics, natural language processing, and other settings in which 
significant portions of the distribution are comprised of domain elements with small probabilities 
(e.g. m)- In its most simple form, the Good-Turing frequency estimation scheme estimates 

the total probability of all domain elements that appear exactly i times in a set of n samples 
as ^ +1 ^ I+1 , where J-j is the total number of species that occur exactly j times in the samples. 
The total probability mass consisting of domain elements that are not seen in the samples—the 
“missing” mass, or, equivalently, the probability that the next sample drawn will be a new domain 
element that has not been seen previously—can be estimated by plugging i = 0 into this formula 
to yield J-\/n , namely the fraction of the samples consisting of domain elements seen exactly once. 

The Good-Turing estimate is especially suited to estimating the total mass of elements that 
appear few times; for more frequently occurring domain elements, this estimate has high variance— 
for example, if Tj+i = 0, as will be the case for most large i, then the estimate is 0. However, for 
frequently occurring domain elements, the empirical distribution will give an accurate estimate of 
their probability mass. There is an extremely long and successful line of work, spanning the past 
60 years, from the computer science, statistics, and information theory communities, proposing 
approaches to “smoothing” the Good-Turing estimates, and combining such smoothed estimates 
with the empirical distribution (e.g. 22 IS I2H E3 12SJ |T5] Sj). 

Our approach—to first recover an estimate of the unlabeled vector of probabilities of the true 
distribution and then assign probabilities to the observed elements informed by this recovered 
vector of probabilities—deviates fundamentally from this previous line of work. This previous work 
attempts to accurately estimate the total probability mass corresponding to the set of domain 
elements observed i times, for each i. Even if one knows these quantities exactly , such knowledge 
does not translate into an optimal learning algorithm, and could result in an t\ error that is a 

x One of the medical benefits of “genome wide association studies” is the compilation of catalogs of rare mutations 
that occur in healthy individuals; these catalogs are being used to rule out genetic causes of illness in patients, and 
help guide doctors to accurate diagnoses (see e.g. |17lll8| l. Understanding how these catalogs will grow as a function 
of the number of genomes sequenced may be an important factor in designing such future datasets. 
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factor of two larger than that of our approach. The following rephrasing of Example [I] from above 
illustrates this point: 

Example 5. Consider n independent draws from a distribution in which 90% of the domain ele¬ 
ments occur with probability 10/n, and the remaining 10% occur with probability 11/n. All variants 
of the Good-Taring frequency estimation scheme would end up, at best, assigning probability 10.1/n 
to most of the domain elements, incurring an t\ error of roughly 0.2. This is because, for elements 
seen roughly 10 times, the scheme would first calculate that the average mass of such elements is 
10.1/n, and then assign this probability to all such elements. Our scheme, on the other hand, would 
realize that approximately 90% of such elements have probability 10 /n, and 10% have probability 
11/n, but then would assign the probability minimizing the expected error — namely, in this case, our 
algorithm would assign the median probability, 10/n, to all such elements, incurring an i\ error of 
approximately 0.1. 

Worst-case vs. Instance Optimal Testing and Learning. Sparked by the seminal work of 
Goldreich, Goldwasser and Ron [20) and that of Batu et al. im there has been a long line of work 
considering distributional property testing, estimation, and learning questions from a worst case 
standpoint—typically parameterized via an upper bound on the support size of the distribution 
from which the samples are drawn (e.g. [HI EH El [23l uni [271 [331 EQl EH [m EU). 

The desire to go beyond this type of worst-case analysis and develop algorithms which provably 
perform better on “easy” distributions has led to two different veins of further work. One vein 
considers different common types of structure that a distribution might possess-such as mono¬ 
tonicity, unimodality, skinny tails, etc., and how such structure can be leveraged to yield improved 
algorithms mum- While this direction is still within the framework of worst-case analysis, the 
emphasis is on developing a more nuanced understanding of why “easy” instances are easy. 

Another vein of very recent work beyond worst-case analysis (of which this paper is an example) 
seeks to develop “instance-optimal” algorithms that are capable of exploiting whatever structure 
is present in the instance. For the problem of identity testing—given the explicit description of p, 
deciding whether a set of samples was drawn from p versus a distribution with i\ distance at least e 
from p —recent work gave an algorithm and an explicit function of p and e that represents the sample 
complexity of this task, for each p [32 ], In a similar spirit, with the dual goals of developing optimal 
algorithms as well as understanding the fundamental limits of when such instance-optimality is not 
possible, Acharya et al. have a line of work from the perspective of competitive analysis mmmm- 
Broadly, this work explores the following question: to what extent can an algorithm perform as 
well as if it knew, a priori, the structure of the problem instance on which it was run? For example, 
the work [2] considers the two-distribution identity testing question: given samples drawn from two 
unknown distributions, p and q, how many samples are required to distinguish the case that p = q 
from ||p —g||i > e? They show that if n P)(? is the number of samples required by an algorithm that 
knows, ahead of time, the unlabeled vector of probabilities of p and q, then the sample complexity 
is bounded by n p \f, and that, in general, a polynomial blowup is necessary—there exists p,q for 
which no algorithm can perform this task using fewer than n p (q samples. 

Relation to [29 . 3J ]. The papers [591 (3Tj were concerned with developing estimators for entropy, 
support size, etc.—properties that depend only on the unlabeled vector of probabilities of a dis¬ 
tribution. The goal in those papers was to give tight worst-case bounds on these estimation tasks 
in terms of a given upper bound on the support size of the distribution in question. In contrast, 
this work considers the question of learning probabilities for each labeled domain element, and 
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considers the more ambitious and practically relevant goal of “instance-optimality”. This present 
paper has two technical components corresponding to the two stages of our algorithm: the first 
component is recovering an approximation to the unlabeled vector of probabilities, and the sec¬ 
ond part leverages the recovered unlabeled vector of probabilities to output a labeled vector. The 
majority of the technical machinery that we employ for the first part is based on algorithms and 
techniques developed in [29, 3TJ, though analyzed here in a more nuanced and general way (a main 
tool from these works is a Chebyshev polynomial earthmover scheme, which was also repurposed 
for a rather different end in m, the main improvement in the analysis is that our results no longer 
require any bound on the support size of the distribution, and the results no longer degrade with 
increasing support size). We are surprised and excited that these techniques, originally developed 
for establishing worst-case optimal algorithms for property estimation can be fruitfully extended 
to yield tight instance-optimal results for such a fundamental and classic learning question. 

1.3 Definitions 

We refer to the unlabeled vector of probabilities of a distribution as the histogram of the distribution. 
This is simply the histogram, in the usual sense of the word, of the vector of probabilities of the 
domain elements. We give a formal definition: 

Definition 1. The histogram of a distribution p, with a finite or countably infinite support, is 
a mapping h p : (0,1] -> NU {0}, where h p (x ) is equal to the number of domain elements that 
occur in distribution p with probability x. Formally, h p (x) = |{a : p(a) = x}\, where p(a) is the 
probability mass that distribution p assigns to domain element a. We will also allow for “generalized 
histograms” in which h p does not necessarily take integral values. 

In analogy with the histogram of a distribution, it will be convenient to have an unlabeled 
representation of the set of samples. We define the fingerprint of a set of samples, which essen¬ 
tially removes all the label-information. We note that in some of the literature, the fingerprint is 
alternately termed the pattern , histogram, histogram of the histogram or collision statistics of the 
samples. 

Definition 2. Given samples X = (x\,.. ., x n ), the associated fingerprint, F = (F\, F 2 , ■ ■ •), is 
the “histogram of the histogram” of the sample. Formally, F is the vector whose i th component, 
Fi, is the number of elements in the domain that occur exactly i times in X. 

The following earthmover metric will be useful for comparing histograms. This metric is sim¬ 
ilar to that leveraged in [29], but allows for discrepancies in sufficiently small probabilities to be 
suppressed. This turns out to be the “right” metric for establishing our learning result, as well 
as our result for the accurate estimation of the expected number of distinct elements that will be 
observed for larger sample sizes (Proposition [l]). In both of these settings, we do not need to worry 
about accurately estimating extremely small probabilities, as long as we can accurately estimate 
the total aggregate probability mass comprised of such elements. 

Definition 3. For two distributions p\,P 2 with respective histograms h\,h 2 , and a real number 
t € [0,1], we define the r-truncated relative earthmover distance between them, R t {pi,P 2 ) ■= 
R r (h \, / 12 ), as the minimum over all schemes of moving the probability mass in the first histogram 
to yield the second histogram, where the cost per unit mass of moving from probability x to probability 
y is | log max(x, r) — log max(y, r) |. 
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The following fact, whose proof is contained in Appendix [Bj relates the r-truncated relative 
earthmover distance between two distributions, p\,P 2 , to an analogous but weaker statement about 
the t\ distance between p\ and a distribution obtained from p 2 by choosing an optimal relabeling 
of the support: 

Fact 1. Given two distributions P\-P 2 satisfying R T (pi,P 2 ) < e, there exists a relabeling n of the 
support of P 2 such that ]>T |max(pi(i), r) — max(p2(7r(*)), r )| < 2e. 

The Poisson distribution will also feature prominently in our algorithms and analysis: 

Definition 4. For A > 0, we define Poi(X) to be the Poisson distribution of parameter A, where 

_A \ 7 

the probability of drawing j <r- Poi( A) equals poi(\,j) = e . 

2 Recovering the histogram 

This section adapts the techniques of [29] [31] to accurately estimate the histogram of the distribution 
in a form that will be useful for Algorithm [2l our ultimate instance-optimal algorithm for learning 
the distribution, presented and analyzed in Section 02 

The first phase of our algorithm, the step in which we recover an accurate approximation of the 
histogram of the distribution from which the samples were drawn, consists of solving an intuitive 
linear program. The variables of the linear program represent the histogram entries h(x i), h(x 2 ),... 
corresponding to a fine discretization of the set of probability values 0 < x\ < X 2 < ■ ■ ■ < 1. The 
constraints of the linear program represent the fact that h corresponds to the histogram of a 
distribution, namely all the probabilities sum to 1, and the histogram entries are non-negative. 
Finally, the objective value of the linear program attempts to ensure that the histogram h output 
by the linear program will have the property that, if the samples had been drawn from a distribution 
with histogram h, the expected number of domain elements observed once, twice, etc., would closely 
match the corresponding actual statistics of the sample. Namely, the objective function tries to 
ensure that the expected fingerprint of the histogram returned by the linear program is close to the 
actual fingerprint of the samples. 

One minor subtlety is that we will only solve this linear program for the portion of the histogram 
corresponding to domain elements that are not seen too many times. For elements seen very 
frequently (at least n a times for some appropriately chosen absolute constant a > 0) their empirical 
probabilities are likely quite accurate, and we simply use these probabilities. A similar approach was 
fruitfully leveraged in [29], 31] with the goal of creating worst-case optimal estimators for entropy, 
and other distributional properties of interest, and a related heuristic was proposed in the 1970’s by 
Efron and Thisted [16], also with the goal of estimating properties of the underlying distribution. 

We state the algorithm and its analysis in terms of two positive constants, £>,C, which can be 
defined arbitrarily provided the following inequalities hold: 0.1 > B > C > § > 0. 
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Algorithm 1. 

Input: Fingerprint T obtained from n-samples. 

Output: Histogram hpp. 

• Define the set X := {—y , -4-, ■ ■ •, n +n }■ 

• For each x £ X, define the associated variable t> x , and solve the LP: 

n B 

Minimize 

i= 1 

Subject to: 

' E xe x ^ + Er>ne +2n c ^ = 1 (total prob. mass = 1) 

• \/x £ X,v x > 0 (histogram entries are non-negative) 

• Let hpp be the histogram formed by setting hLp(xi) = v Xi for all i, where (v x ) is the solution 
to the linear program, and then for each integer i > nr + 2n , incrementing hLp(-) by T % . 


Ti¬ 


ll 

x&X 


poi(nx, i) ■ v a 


The following theorem quantifies the performance of the above algorithm: 


Theorem 2. There exists an absolute constant c such that for sufficiently large n and any w £ 
[l,logn], given n independent draws from a distribution p with histogram h, with probability 1 — 
e~ n the generalized histogram hpp returned by Algorithm]]} satisfies 


/? m 

n log n 


(h, hLp) < —j=- 
1 Jw 


This theorem is a stronger and more refined version of the results in (29] , in that these results no 
longer require any bound on the support size of the distribution, and the results no longer degrade 
with increasing support size. Instead, we express our results in terms of a lower bound, r, on the 
probabilities that we are concerned with accurately reconstructing. We provide the proof of the 
theorem in a self-contained fashion in Appendix [Dj 

One interpretation of Theorem [2] is that Algorithm |T] when run on n independent draws from a 
distribution, will accurately reconstruct the histogram, in the relative earthmover sense, all the way 
down to probability nl( ^ gn (significantly below the 1/n threshold where the empirical distribution 
becomes ineffective). One corollary of independent interest is that this earthmover bound implies 
that we can accurately extrapolate the number of unique elements that will be seen if we run a 
new, larger experiment, of size up to nlogn. Given a histogram h, for each element of probability 
x , the probability that it will be seen (at least once) in a sample of size k equals 1 — (1 — x) k \ 
thus, the expected number of unique elements seen in a sample of size k for a distribution with 
histogram h equals XE/i(x)^o(l ~ (1 ~ x ) k ) ' h{x). The following lemma, whose proof is given in 
Appendix O shows that this quantity is Lipschitz continuous with respect to truncated relative 
earthmover distance. 


Lemma 1. Given two (possibly generalized) histograms g, h, a number of samples k, and a threshold 








TG (0,1] , 


< (0.3 (k - 1) + 1 )Rr(g, h) + 


E (i-(i-E)-sW- E (i-(i-*)*)■ aw 

x:g(x)^0 x:h(x)^ 0 

The above lemma together with Theorem [2] yields Proposition [Q which is tight, in the sense 
that one cannot expect meaningful extrapolation beyond sample sizes of nlogn, as shown by the 
lower bounds in [29] H 

Towards our goal of devising an optimal learning algorithm, the following corollary of Theorem [2] 
formalizes the sense that the quality of the histogram output by Algorithm [T] will be sufficient to 
achieve our optimal learning result, provided that the second phase of our algorithm, described in 
Section [31 is able to successful label the histogram. 


Corollary 1. There exists an algorithm such that, for any function f{n) = w n (l) that goes to 
infinity as n gets large (e.g. f(n) = log log n), there is a function o n { 1) of n that goes to zero as 
n gets large, such that given n samples drawn independently from any distribution p, the algorithm 
outputs an unlabeled vector, q, such that, with probability 1 — e -T?n(1) , there exists a labeling ir(q) of 
the vector q such that 



fin) \ 
n log n J 


— max 



f{n) \ 
n log n ) 


o n (l), 


where p(x) denotes the true probability of domain element x in distribution p. 

This corollary is not immediate: the histogram returned by the algorithm might be non-integral, 
however in Appendix [E] we provide a simple algorithm that rounds a generalized histogram to an 
(integral) histogram, while changing it very little in relative earthmover distance This 

rounding, together with Fact [H obtains this corollary. 

The utility of the above corollary lies in the following observation: for any function g(n) = 
o(l/n), the domain elements x that both occur in the n samples and have true probability p(x) < 
g{n), can account for at most o(l) probability mass, in aggregate. In other words, while the true 
distribution might have a constant amount, c, of probability mass consisting of domain elements 
that occur with probability o(l/n), we would observe at most a o(l) fraction of such domain 
elements in the n samples. Hence, even an optimal scheme that knows the true probabilities would 
be unable to achieve an t\ error less than c — o(l) because it does not know the labels of the 
elements that have not been observed, and we could also hope to achieve an t\ error of roughly c. 


3 Disambiguating the histogram 

In this section we present our instance-optimal algorithm for learning a distribution from n samples, 
making use of Algorithm |T] of Section [2] to first accurately infer the histogram of the distribution 

2 Namely, for some constant c, there exist two families of distribution, X>i and £>2 such that the distributions in 
X>i are close to uniform distributions on cn log n elements, and the distributions in T >2 are close to uniform distribu¬ 
tions over 2cnlogn elements, yet a randomly selected distribution in T> 1 is (with constant probability) information 
theoretically indistinguishable from a randomly selected distribution in T >2 if one is given only n samples drawn from 
the distributions. Of course, given 2cnlogn samples, the number of distinct elements observed will likely be either 
~ (1 — -fi)cn log n « 0.9cnlogn or « (1 — i)2cnlogn « 1.3cnlogn according to the two cases. 
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(in the sense of Corollary [TJ). As a motivating intuition for the second phase of our algorithm— 
the phase in which we assign probabilities to the observed elements—consider the behavior of an 
optimal algorithm that not only knows the true histogram h of the distribution, but also knows for 
each positive integer j the entire multiset of probabilities of elements that appear exactly j times 
in the n samples. Since the algorithm has no basis to distinguish between the different elements 
that each appear j times in the samples, the algorithm may as well assign a single probability rrij 
to all the items that appear j times in the samples. The optimal mj in this setting is easily seen 
to be the median of the multiset of probabilities of items appearing j times, as the median is the 
estimate that minimizes the total (£ 1 ) error of the probabilities. 

Our algorithm aims to emulate this idealized optimal algorithm. Of course, we must do this 
using only an estimate of the histogram, and computing medians based on the likelihoods that 
elements of probability x will be seen j times in the sample, as opposed to actual knowledge of 
the multiset of probabilities of the elements observed j times (which was an unreasonably strong 
assumption, that we made in the previous paragraph because it let us argue about the behavior of 
the optimal algorithm in that case). 

Because our algorithm needs to work in terms of a histogram estimate u , bounded only by the 
guarantees of Corollary |TJ we add an additional “regularization” step that was not needed in the 
idealized medians setting described above. We “fatten” the histogram u to a new histogram u by 
adding a small amount of probability mass across the range ^ log 2 n], which acts to mollify the 
effect on the medians of any small errors in the histogram estimate. 

Given this “fattened” approximate histogram, we then apply the “medians” intuition to com¬ 
puting, for each integer j. an appropriate probability with which to label those elements occurring j 
times in the sample. These estimates are computed via the following thought experiment: imagining 
u to be the true histogram, if we take n samples from the corresponding distribution, in expectation , 
what is the median of the (multiset) of probabilities of those elements seen exactly j times in the 
sample? We denote this “expected median” by my j, n . and our algorithm assigns this probability to 
each element seen j times in the sample, for j < log 2 n, and assigns the empirical probability - for 
larger j. We formalize this process with the following definition for “Poisson-weighted medians”: 

Definition 5. Given a histogram h, let Sh be the multiset of probabilities of domain elements—that 
is, for each probability x for which h(x ) is some positive integer i, add i copies of x to S. Given a 
number of samples n, and an index j, consider weighting each element x € Sh by poi{nx,j). Define 
m h,j,n t° be the median of this weighted multiset,. 

Explicitly, the median of a weighted set of real numbers is a number m, such that at most half 
the weight lies on numbers greater than m, and at most half lies on numbers less than rn. Taking 
advantage of the medians defined above, our learning algorithm follows: 
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Algorithm 2. 

Input: n samples from a distribution h. 

Output: An assignment of a probability to each nonzero entry of h. 

• Run Algorithm [T] to return a histogram u. 

• Modify u to create u by, for each j < log 2 n adding lo " 4 - elements of probability ^ and 
removing corresponding mass arbitrarily from the rest of the distribution. 

• Then to each fingerprint entry j < log 2 n, assign those domain elements probability mu,j,n, 
(as defined in Definition [5]) and to each higher fingerprint entry j > log 2 n assign those domain 
elements their empirical probability A 


Theorem |T] There is a function err{n ) that goes to zero as n gets large, such that Algorithm 
when given as input n independent draws from any distribution p of discrete support, outputs a 
labeled vector q, such that 

E[\\p — (?||i] < opt(p,n) + err(n), 

where opt(p , n) is the minimum expected error that any algorithm could achieve on the following 
learning task: given p, and given n samples drawn independently from a distribution that is identical 
to p up to an arbitrary relabeling of the domain elements, learn the distribution. 

The core of the proof of Theorem [T] relies on constructing an estimate, devj^ n {A,mB,j, n )i that 
captures the expected contribution to the i\ error due to elements that occur exactly j times, given 
that the true distribution we are trying to reconstruct has histogram A, and our reconstruction is 
based on the medians mB,j, n derived from a (possibly different) histogram B. The proof then has 
two main components. First we show that devj_ n {h, mh,j, n ) approximately captures the performance 
of the optimal algorithm with very high probability, namely that using the true histogram h to 
choose medians mh } j, n lets us estimate the performance of the best possible algorithm. This step is 
slightly subtle, and implies that an algorithm that knows h can glean at most o(l) added benefit by 
computing the probability assigned to an element that occurs j times using a function that depends 
on j and h and the entire set of samples, rather than just j and h. 

Next, we show that the clean functional form of this “median” estimate implies that dev(-,-) 
varies slowly with respect to changes in the second histogram (used to choose the median in the 
second term), and thus that with only negligible performance loss we may reconstruct distributions 
using medians derived from an estimate u of the true histogram, thus allowing us to analyze the 
actual performance of Algorithm [2j 

Beyond these two core steps, the analysis of Algorithm [2] is somewhat delicate—because our al¬ 
gorithm is instance-optimal to o(l) error, it must reuse samples both for the Algorithm Q] histogram 
reconstruction and for the final labeling step, and we must carefully separate the probabilistic por¬ 
tion of the analysis via a clean set of assumptions which 1) will hold with near certainty over the 
sampling process, and 2) are sufficient to guarantee the performance of both stages of our algorithm. 
The complete proof is contained in Appendix lAl 
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A Proof of Theorem [T] 

In this section we give a self-contained proof of the correctness of Algorithm [2], establishing Theo¬ 
rem in 

A.l Separating out the probabilistic portion of the analysis 

Our analysis is somewhat delicate because we reuse the same samples both to estimate the histogram 
h, and then to label the domain elements given an approximate histogram. For this reason, we will 
very carefully separate out the probabilistic portion of the sampling process, identifying a list of 
convenient properties which happen with very high probability in the sampling process, and then 
deterministically analyze the case when these properties hold, which we will refer to as a “faithful” 
set S of samples from the distribution. (Appendix [D] uses this same analysis technique, though 
with a different notion of “faithful”, appropriate for the different desiderata of that appendix.) 
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We first describe a simple discretization of histograms h, dividing the domain into buckets which 
will simplify further analysis, and is a crucial component of the definition of “faithful”. 

Definition 6. Given a histogram h, and a number of samples n, define the kth bucket ofh to consist 
of those histogram entries with probabilities in the half-open interval ( nlo g^ n ; n f 0 g^ n ]- Letting hk 
be h restricted to its kth bucket, define B po fij,k) = J2x-h k (x)^oh( x )p°i( nx ij) to be the expected 
number of elements from bucket k that are seen exactly j times, if Poi(n) samples are taken. Given 
a set of samples S, let Bs(j, k ) be the number of elements in bucket k of h that are seen exactly j 
times in the samples S, where in both cases h and n are implicit in the notation. 

Given this notion of “buckets”, we define faithful to mean 1) each domain element is seen 
roughly the number of times we would expect to see it, and 2) for each pair ( j,k ), the number 
of domain elements from bucket k that are seen exactly j times is very close to its expectation 
(where we compute expectations under a Poisson distribution of samples, because “Poissonization” 
will simplify subsequent analysis). The first condition of “faithful” gives weak control on which 
fingerprint entry each domain element will contribute to, while the second condition gives much 
stronger control over the aggregate contribution to fingerprint entries by all domain elements within 
a certain probability “bucket”. 

Definition 7. Given a histogram h and a number of samples n, a set of n samples, S, is called 
faithful if: 

1. Each item of probability x appears in the samples a number of times j satisfying \nx — j\ < 
maxjlog 1 ' 5 n, \Jnx log 1 " 5 n}, and 

2. For each j < log 2 n and k, we have \B po fij , k ) — Bs(j, k)\ < n 0 ' 6 . 

This notion of “faithful” holds with near certainty, as shown in the following lemma, allowing 
us to assume (when specified) in the results in the rest of this section that our learning algorithm 
receives a faithful set of samples. 

Lemma 2. For any histogram h and number of samples n, with probability 1 — , a set of n 

samples drawn from h will be faithful. 

Proof. Since the number of times an item of probability x shows up in n samples is the binomial 
distribution Bin(n, x), the first condition of “faithful”—essentially that this random variable will be 
within log 3 / 4 n standard deviations of its mean— follows with probability 1 — n - ^ 1 ) from standard 
Chernoff/Hoeffding bounds. 

For the second condition, since Poi(n ) has probability 0(1 /y/n) of equaling n, we consider the 
related process where Poi(n ) samples are drawn. The number of times each domain element x is 
seen is now distributed as Poi(nx), independent of each other domain element. Thus the number of 
elements from bucket k seen exactly j times is the sum of independent Bernoulli random variables, 
one for each domain element in bucket k. The expected number of such elements is B poi (j, k) by 
definition. Since B po fij,k ) < n by definition, we have that the variance of this random variable is 
also at most n, and thus Chernoff/Hoeffding bounds imply that the probability that it deviates from 
its expectation by more than n 0 ' 6 is at most exp(— n ai ). Thus the probability of such a deviation 
is at most a @(y/n) factor higher when taking exactly n samples than when taking Poifn) samples; 
taking a union bound over all j and k yields the desired result. □ 
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A.2 An estimate of the optimal error 

We now introduce the key definition of dev(-,-), which underpins our analysis of the error of 
estimation algorithms. The definition of dev(-,-) captures the following process: Suppose we have 
a probability value rrij . and will assign this probability value to every domain element that occurs 
exactly j times in the samples. We estimate the expected error of this reconstruction, in terms of 
the probability that each domain element shows up exactly j times. While the below definition, 
stated in terms of a Poisson process, is neither clearly related to the optimal error opt(h, n), nor the 
actual error of any specific algorithm, it has particularly clean properties which will help us show 
that it can be related to both opt(h, n ) (in this subsection) as well as the expected error achieved 
by Algorithm [2] (shown in Section IA.3D . 

Definition 8. Given a histogram h, a real number m, a number of samples n, and a nonnegative 
integer j, define devj tn (h,m ) = Yl x -h(x)^o l x — m\h(x)poi(nx,j). 

Intuitively, devj tTl (h. rn) describes the expectation—over taking Poi(n) samples from h —of the 
sum of the deviations between m and each probability x of a element seen j times among the 
samples. Namely, devj <n (h,m ) describes to what degree rn is a good probability to which we can 
ascribe all domain elements seen j times, among ~ n samples from h. 

This definition provides crucial motivation for how Definition [5] sets the medians rri^ j n used in 
Algorithm^ since rri}- L ] ri is the value of m that minimizes the previous definition, devj tn (h , m ), since 
both are defined via the same Poisson weights poi(nx,j). (The median of a- possibly weighted— 
set of numbers is the location m that minimizes the total—possibly weighted—distance from the 
set to m.) 

We now show the key result of this section, that the definition of “faithful” induces precise 
guarantees on the spread of probabilities of those elements seen j times. Subsequent lemmas will 
relate this to the performance of both the optimal algorithm and to our own Algorithm [2j 

Lemma 3. Given a histogram h, let S be the multiset of probabilities of a faithfid set of samples 
of size n. For each index j < log 2 n, consider those domain elements that occur exactly j times 
in the samples and let Sj be the multiset of probabilities of those domain elements. Let aj be 
the sum over Sj of each element’s distance from the median (counting multiplicity) of Sj. Then 
Ej<io g 2 n Wj - dev jin (h,m h j jn )\ = 0(log“ 2 n). 

Proof. Recall that a computes the total distance of the (unweighted) multiset Sj from its median, 
while devj in (h,mh,j t n ) is an analogous (weighted) quantity for the true histogram, with each entry 
x having multiplicity h(x ) and weight poi(nx,j). In the first case, sampling means that each 
element of probability x either shows up exactly j times (with some binomial probability) and is 
counted with weight 1, or does not show up j times and is not counted; in the second case, instead of 
sampling, each entry x from the histogram is counted with weight poi(nx,j) < 1, capturing roughly 
the average effect of sampling (except with Poisson instead of binomial weight). By the definition 
of “faithful”, the total weight coming from each bucket k in both cases is within n 0 ' 6 of each other 
(since j < log 2 n). We consider only buckets k < 2 log 2 ' 2 n, corresponding to probabilities less than 
— log 2 n, since the first condition of “faithful” means that no higher probability elements will be 
seen j < log 2 n times. 

Consider transforming one weighted multiset into the other (where elements of Sj are interpreted 
as having weight 1 each), maintaining a bound on how much the total distance from the median 
changes. We make crucial use of the fact that the “total distance to the median” is robust to small 
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changes in the weighted multiset, since the median is the location that minimizes this total distance. 
Moving a weight by a distance of (3 can increase the total (weighted) distance to the median by at 
most a ■ (3 since this is how much the total weighted distance to the old median changes, and the 
new median must be at least as good; conversely, such a move cannot decrease the total distance 
by more than a ■ f3 as the inverse move would violate the previous bound. Adding a weight to the 
distribution at distance /3 from the current median similarly cannot decrease the total distance, 
but also cannot increase the total distance by more than a ■ (3, with the corresponding statements 
holding for removing a weight. 

Thus, transforming all the Sj into their Poissonized analogs requires two types of transforma¬ 
tions: 1) moving up to n samples within their buckets; 2) adding or removing up to n 0 ' 6 weight from 
buckets for various combinations of j and k. Since buckets have width l/(nlog 2 n), transforma¬ 
tions of the first type change the total distance to the median by at most log -2 n; since j < log 2 n 
and all buckets above probability ^ log 2 n are empty, transformations of the second type change 
the total distance by at most the product of the weight adjustment n a6 , the number of j,k pairs 
2 log 3 ' 2 n, and size of the probability range under consideration which is ^ log 2 n, yielding a bound 
of log 4 ' 2 n. Thus in total the change is 0(log -2 n) = o(l) as desired. □ 

The above lemma essentially shows that devj^ n (h,mhj^ n ) captures how well we could hypothet¬ 
ically estimate the probabilities of all the domain elements seen j times, under the unrealistically 
optimistic assumption that we know the (unlabeled) multiset of probabilities of elements seen j 
times and estimate all these probabilities optimally by their median. Before showing how our 
algorithm can perform almost this well based on only the samples, we first formalize this reasoning. 

Definition 9. We call a distribution learner “simple” if all the domain elements seen exactly j 
times in the samples get assigned the same probability. 

Given n samples from a distribution p, with being those domain elements that occurred 
exactly j times in the sample, we note that the probability of obtaining these samples is invariant to 
any permutation of P(j). Thus if a hypothetical learner L assigns different probabilities to different 
elements seen j times in the sample, then its average performance over a random permutation 
of the domain elements can only improve if we simplify L, by having it instead assign to all the 
elements seen j times the median of the multiset that it was originally assigning. 

For this reason, when we are discussing an optimal distribution learner, we will henceforth 
assume it is simple. 

Lemma 4. Given a histogram h, let S be the multiset of probabilities of a faithful set of samples 
of size n. Given an index j < log 2 n, consider those domain elements that occur exactly j times 
in the sample; let Sj be the multiset of probabilities of those domain elements. Let a j be the sum 
over Sj of each element’s distance from the median of Sj (counting multiplicity). Then any simple 
learner, when given the sample, must have error at least Oj on the domain elements that appear j 
times in the sample. 

Proof. The median of Sj is the best possible estimate any simple learner can yield—even given the 
true distribution—so the error of this estimate bounds the performance of a simple learner. □ 

Combining this with Lemma [3] immediately yields: 
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Corollary 2. For any distribution h, the total error of any simple learning algorithm, given n 
faithful samples from h, is at least — 0(log _2 n). Further, for any 

algorithm—simple or not—if we average its performance over all relabelings of the domain of h and 
the corresponding relabeled samples, it will have expected error bounded by the same expression. 

A.3 Our error estimate is Lipschitz with respect to mis-estimating the distri¬ 
bution 

We now relate the error bound of Corollary [2] to the performance of our algorithm, via two steps. 
The bound in the corollary is in terms of mh,j,ni the medians computed in terms of the true 
histogram h which is unknown to the algorithm; instead the algorithm works with an estimate u 
of the true histogram. The next lemma shows that estimating in terms of u is almost as good as 
using h. 

Fact 2. For any distribution h, index j > 1, and real parameter t > 1, weighting each domain 
element x by poi(nx,j), the total weight on domain elements that are at least t standard deviations 
away from £— namely, for which \nx — j\ > ty/j is at most n ■ exp(—Q(t)). 

Lemma 5. Given a number of samples n, a histogram h and a second histogram u that is 1) close to 
h in the sense of Corollary [If in that there exists distributions p, q corresponding to h,u respectively 
for which JT | ma x(p(i), -t log -0 ’ 25 n) — ma x(q(i), ^ log -0 ' 25 n)\ < log -0 ' 25 n, and 2) the histogram u 
is “ fattened” in the sense that for each j < log 2 n there are at least . lo ^.a — elements of probability 

n- Tken D j<\og 2 n deV j,n( h ’ m u,j,n ) < o(l) + Ej<log 2 n dev j,n(h, ■ 

Since for each j, as noted earlier, mh,j, n is the quantity which minimizes devj tU (h, m), each term 
devj.rf h, riuj^n) on the left hand side is greater than or equal to the corresponding devj_ n (h, n^n) 
on the right hand side, so the lemma implies that the left and right hand sides of the expression in 
the lemma, beyond having related sums, are in fact term-by-term close to each other. 

The proof relies on first comparing and muj to A and then showing that devj(h,m) is 
Lipshitz with respect to changes in h of the type described by the guarantees of Corollary [lj 

Proof of Lemma El We drop the n” subscripts here for notational convenience. 

Recall that the quantities rrih.j and are medians computed after weighting by a Poisson 
function centered at j, and thus we would expect these medians to be close to A We first show 
that the “fattening” condition makes well-behaved, and then show, given this, that the lemma 
works both in the case that mhj is far from j-. and then for the case where they are close. 

By condition 2 of the lemma, the “fattening” assumption, for any index j < log 2 n, we have 
h(x)poi(nx , j) = l/log°^ n. Thus, by FactEl the median must satisfy \n-m.uj—j\ < 
y/j log^ 1 - 1 n, since the fraction of the Poisson-weighted distribution that is at locations more than 
y t yfj log 0 ^ n distance from - i s (much) less than 1/2. 

Given the above bound on m S j-, we now turn to mh,j- Consider the case |n ■ mhj — j\ > 
y/jlog 01 n. By Fact [2l weighting each domain element x by poi(nx,j), the total weight on the 
far side of the median m^j from is at most n ■ exp(— ^(log 0 ' 1 n)). Since (by definition of 
“median”) half the weight is on each side of the median, the total weight Yh x -h(x)^o h(x)p°i(nx, j) 
must also be bounded by n ■ exp{—Ll{ log 01 n)). Recall the definition of the left hand side of 
the inequality of the lemma, devj{h,niuj ) = Yl x -h(x)^o \ x ~ rriu l j\h(x)poi(nx, j). Thus for the 
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portion of this sum where x < y log 2 n, since from the previous paragraph is also bounded by 
^log 2 n for large enough n, we can bound Yl x <?-iog 2 n -h.( x )^o \ x ~ m u,j\h{x)poi(nx,j) by the product 
n ■ exp(—td(log 01 n)) ■ ^log 2 n = exp(— n(log ai n)). For those x > |log 2 ?r, since j < ^ log 2 n, 
we have the tail bounds poi(nx,j ) = , implying the total for such x is also bounded by 

exp(—Q( log 0 ' 1 n)), which is our final bound for this case—summing these bounds over all j < log 2 n 
yields the desired bound X^<iog 2 n devj(h, m u,j) < o(l), where the sum is over those j for which 
this case applies, | n ■ rn.hj — j\ > yfj log 0 ' 1 n. 

Thus it remains to prove the claim when both and rriuj are close to /■. To analyze this 
case, we show that devj(h,m) is Lipschitz with respect to the closeness in h and u guaranteed by 
condition 1 of the lemma, provided \n ■ m — j\ < y/j log 01 n. The guarantee on h and u means that 
one can transform one distribution into the other by two kinds of transformations: 1) changing 
the distributions by log -0 ' 25 n in the i\ sense, and 2) arbitrary mass-preserving transformation of 
elements of probability less than ^ log^ 0 ' 25 n. We thus bound the change in devj(h, rn) under both 
types of transformations. 

To analyze t\ modifications, consider an arbitrary probability x, and consider the derivative 
of devj(h,m ) as we take an element of probability x and change x. Recalling the definition 
devj(h,m ) = Y/ x -h(x)^o I x — m\h(x)poi(nx,j), we see that this derivative equals -^\x — m\poi(nx, j), 
which is bounded (by the product rule and triangle inequality) as poi(nx,j) + \x — m\ ■4-poi(nx,j), 
where /^poi(nx, j) = n ■ poi(nx,j — 1) ■ (1 — ^). Rewriting m as mj to indicate its depen¬ 
dence on j, we want to bound the sum of this derivative over j < log 2 n, since the exact depen¬ 
dence for each individual j is much harder to talk about than the overall dependence. We have 
^2jPoi(nx,j) + \x — rrij\n-poi(nx,j — 1) • (1 — ^), where Y/jpoi(nx,j) < 1. To bound the remaining 
part of the sum, we first consider the case x < in which case we bound \x — rn 3 \ < ^(1 + j + 
yfj log 01 n) and (1 — < 1, thus yielding the bound X^'>i \ x ~ m j\ n ' poi(nx,j — 1) • (1 — < 

Ej>o ( 2 +J + Vd + 1 log 0 ' 1 n)poi(nx, j ) < Yj>o( 2 +d+ VJ + Hog 0 ' 1 n)/j\ = 0(log ai n). For x > 
since poi(nx,j — 1) decays exponentially fast for j more than yfnx away from nx, we can bound 
this sum as being on the order of yfnx times its maximum value when j is in this range. In this 
range we have \x — nij\ < \x — — m,j\ = ^0(yfnx log 0 ' 1 n ), and poi(nx,j — 1) = 0(-j=), 

and (1 — ?j-) = 0(1/yfnx), yielding a total bound of 0(yfnxyfnx^=^^ log 0 ' 1 n) = ©(log 0 ' 1 n) as 
in the previous case. Thus we conclude that the sum over all j of the amount devj(h,m) changes 
with respect to t\ changes in h is ©(log 0 ' 1 n). 

We next bound the total change to devj(h,m) induced by the second type of modification, 
arbitrary mass-preserving transformations of elements of probability x < ^ log -0 ' 25 n. For j = 1, 
we bound the components of devj(h,m ) = Y x -h(x)^o \ x ~ m\h(x)poi(nx, j) by bounding the two 
terms in the product: \x — m\ G [rn — f log -0 ' 25 n, m + ^ log -0 ' 25 n], and poi(nx, 1) = nx ■ e~ nx G 
[nx( 1 — log -0 ' 25 n) 2 , nx]. Thus for m either m/ ( , or mu, since by the assumption of this case m < 
^(1Tlog 0 ' 1 n), from the bounds above, the contribution to devj(h, m) from those x < ^ log^ 0 ' 25 n is 
within o(l) of mn times the total mass in the distribution below /- log -0 ' 25 n, showing that arbitrary 
modifications of the second type modify dev\(h,m ) by o(l). 

Analyzing the remaining j > 2 terms, omitting the \x — m\ multiplier for the moment, we have 
that SkI io g -°- 25 n-h(x)^o h(x)poi(nx, j) < n( log -0 ' 25 n)- 7-1 . Because of the bound that rri),, rriu are 

each within ^y/Jlog °' 1 n of we have that \x — m\ < ^(log^°' 25 ?r + j + v^Jlog 0 ' 1 ^,). Thus the 
change to devjih, rn) from changes of the second type, summed over all j > 2, is bounded by the 
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sum ^ j>2 (log-°- 25 n )- 7_1 (log -0 ' 25 n + j + yfj log 0 ' 1 n) = o(l), as desired. 

Putting the pieces together, the closeness of h and u implies by the above Lipschitz argument 
that changing the distribution between h and u. under the fixed median niuj does not increase 
dev(-,-) too much: X^ciog 2 n devj(h, m u ,j) < o(l) + Z]j<i og 2 n devj(u, m^j)- Further, since m n j 
minimizes this last expression, the right hand side can only increase if we replace devj(u,rriuj ) by 
devj(u,mhj ) in this last inequality. Finally, a second application of the same Lipschitz property 
implies Ylj<\og 2 n devj{u , rrihj) < o(l) + X^<iog 2 n devj(h, m^j). Combining these three inequalities 
yields the bound of the lemma, X^dog 2 n de y j(h, rriuj) < o(l) + X^j<iog 2 n devj(h, m h,j), as desired. 

□ 

The following lemma characterizes the effect of “fattening” in the second step of Algorithm [2 
showing that this slight modification to the histogram keeps the resulting medians small enough 
that me may apply the following Lemma [71 

Lemma 6. For sufficiently large n, given a fattened distribution fi, for any j < log 2 n, the median 
m p,,j,n is at most 2 log 2 n. 

Proof. Recall that mp^ n is defined as the median of the multiset of probabilities of u after each 
probability x has been weighted by poi(xn, j). For x > j log 2 n and j < log 2 n, these weights will 
each be small by Poisson tail bounds; and because of the fattening, the elements added at 

probability - will contribute inverse poly logarithmic weight. Since the median must have at most 
half the weight to its left, the median cannot be as large as our bound y log 2 n, as desired. □ 

Given the above bound on the size of medians for small j, the following lemma shows that 
our dev(-,-) estimates accurately capture the performance of these medians on any faithful set of 
samples. 

Lemma 7. Given a histogram h, a number of samples n, and for each fingerprint entry j < log 2 n 
a probability mj < ^ log 2 n to which we attribute each domain element that shows up j times in the 
sample, then for any faithful set of samples from h, the total error made for all j < log 2 n is within 
°( 1 ) o/Ej<log 2 n det, J>(/l,mj). 

Proof. Recalling the “buckets” from Definition [6l consider for arbitrary integer k, those elements 
of h in bucket k, which we denote hp —namely, those probabilities of h lying in the interval 
( , - . f' +1 1. where by the first condition of “faithful”, none of these probabilities are above 

^ log 2 n for large enough n. Further, let S-jp be the multiset of probabilities of those domain ele¬ 
ments from bucket k of h that each get seen exactly j times in the sample. The total error of our 
estimate rrij on bucket k is thus YlxeSj k \ m i ~ x l> w hich since buckets have width l/(nlog 2 n), is 
within |5j i fc|/(nlog 2 n) of \Sj t k\ ■ |rrij — A:/(nlog 2 n)|, where we have approximated each x by the left 
endpoint of the bucket containing x. By the second condition of “faithful”, Sjj. is within n 0 6 of 
its expectation, B po i(j, k), and since by assumption mj < log 2 n, we have that our previous error 
bound IS^fcl • |mj — A:/(nlog 2 n)| is within log 2 n of B po i{j,k) ■ \mj — k/(nlog 2 n)\. We rewrite 
this final expression via the definition of B po i as Yl x -h k (x )^o \ m ~ k/(nlog 2 n)\h{x)poi{nx, j). We 
compare this final expression to the portion of the deviation devj^ n (h, mj) that comes from bucket k, 
namely E i:flt (i)/o \mj-x\h{x)poi(nx, j), where since \mj-x\h(x)poi(nx,j) = B poi (j , k) 

and x is within l/(nlog 2 n) of k/(nlog 2 n), the difference between them is clearly bounded by 
B po i(j,k )/(nlog 2 n). Using the triangle inequality to add up the three error terms we have accrued 
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yields that our estimate for the i\ error we make for elements seen j times from bucket k is accurate 
to within 

2 

\Sj,k\/{n\og 2 n) + -^jlog 2 n +B poi (j,k)/(n log 2 n). 

We sum this error bound over all 2 log 2 ' 2 n buckets k and all indices j < log 2 n. The middle 
term log 2 n clearly sums up to o(l) over all j,k pairs. Further, since Sj^ is within n 0 ' 6 of 
Bpoi(jik) by the definition of faithful, the sum of the first term is within o(l) of the sum of the 
third term and it remains only to analyze the third term involving B po i(j,k). From its definition, 
Yj k k) is the expected number of distinct items seen, when making Poi{n) draws from the 

distribution, throwing out those elements which violate the j and k constraints; hence this sum 
over all j, k pairs is at most n, bounding the total error of our u dev v estimates by 0(1/log 2 n), as 
desired. □ 

A.4 Proof of Theorem [I] 

We now assemble the pieces and prove Theorem [lj 

Proof of Theorem QJ Consider the output of Algorithm |T] as run in the first step of Algorithm [2j 
Corollary |T] outlines two cases: with o(l) probability the closeness property outlined in the proposi¬ 
tion fails to hold, and in this case, Algorithm [2] may output a distribution up to t\ distance 2 from 
the true distribution; because this is a low-probability event, this contributes 2 • o(l) = o(l) to the 
expected error. Otherwise, u is close to h, and the fattened version u is similarly close, which lets 
us apply Lemma [5] to conclude that Y^j^og 2 n dev j,n{h, m u,j,n) < o(l) + Yj<\og 2 ndeVj, n (h,m h j tn ). 
Corollary [2] says that ^T <log 2 n devj t71 (h, essentially lowerbounds the optimal error opt(h, n), 

which we combine with the previous bound to yield 

^ devj tn (h, mu,j,n) < opt(h, n) + o(l). 

j<log 2 n 

Lemma [2] guarantees that the samples will be faithful except with o(l) probability, which, as 
above, means that even if these unfaithful cases contribute the maximum possible distance 2 to 
the i\ error, the expected contribution from these cases is still o(l), and thus we will assume 
a faithful set of samples below. Lemmas [6] and [7] imply that for any faithful sample, the error 
made by Algorithm [2] on attributing those elements seen fewer than log 2 n times is within o(l) of 
Yj^nd^Ah,^), and hence at most o(l) worse than opt(h,n). 

Condition 1 of the definition of faithful (Definition [7]) implies that all of the elements seen at 
least log 2 n times originally had probability at least ^(log 2 n — log 175 n) and that the relative error 
between the number of times each of these elements is seen and its expectation is thus at most 
log~ 1//4 ?r. Thus using the empirical estimate on those elements appearing at least log 2 n times—as 
Algorithm [2] does—contributes 0(log _1 ^ 4 n) total error on these elements. Thus all the sources of 
error add up to at most o(l) worse than opt(h,n ) in expectation, yielding the theorem. □ 

B Proof of Fact |T] 

For convenience, we restate Fact [Tj 
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Fact \ 1 \ Given two distributions pi,P 2 satisfying R T (pi ,p 2 ) < e, there exists a relabeling ir of the 
support of P 2 such that 

^2 |max(pi(i),r) - max(p 2 (7r(*)),r)| < 2e. 

i 


Proof of Fact [Jl We relate relative earthmover distance to the minimum L\ distance between re- 
labled histograms, with a proof that extends to the case where both distances are defined above a 
cutoff threshold r. The main idea is to point out that “minimum rearranged” L\ distance can be 
expressed in a very similar form to earthmover distance. Given two histograms hi, /i 2 , the minimum 
L\ distance between any labelings of h\ and /i 2 is clearly the L\ distance between the labelings 
where we match up elements of the two histograms in sorted order. Further, this is seen to equal the 
(regular, not relative) earthmover distance between the histograms hi and /i 2 , where we consider 
there to be hi(x) “histogram mass” at each location x (instead of hi{x) ■ x “probability mass” as 
we did for relative earthmover distance), and place extra histogram entries at 0 as needed so the 
two histograms have the same total mass. 

Given this correspondence, consider an optimal relative earthmoving scheme between hi and 
/i 2 , and in particular, consider an arbitrary component of this scheme, where some probability mass 
a gets moved from some location x in one of the distributions to some location y in the other, at 
cost a log , and suppose without loss of generality that x > y. 

We now reinterpret this move in the Li sense, translating from moving probability mass to 
moving histogram mass. In the non-relative earthmover problem, a probability mass at location x 
corresponds to ^ “histogram mass” at x, which we then move to y at cost (max(r, r) — max(y, r))^; 
however, to simulate the relative earthmover scheme, we need the full ^ mass to appear at y, so 
we move the remaining ^ f mass up from 0, at cost (^ — ^)(max(y,r) — r). 

To relate these 3 costs (the original relative earthmover cost, and the two components of the 
non-relative histogram earthmover cost), we note that if both x and y are less than or equal to r 
then all 3 costs are 0. Otherwise, if x,y > r then the first component of the histogram cost equals 
(1 — |)a and the second is bounded by this, as (^ — ^)(max(y,r) — T ) < (§ — §)y = (1 — Da- 
Further, for the case under consideration where t < y < x, we have (1 — < alog|, which 

equals the relative earthmover cost. Thus the histogram cost in this case is at most twice the 
relative earthmover cost. 

In the remaining case, y < r < x, and the second component of the histogram cost equals 0 
because max(y, r) —r = 0. The first component simplifies as (max(x, r) — max(y, t))j = (x — r)^ = 
(1 — D a < a log D where this last expression is the relative earthmover cost. Thus in all cases, the 
histogram cost is at most twice the relative earthmoving cost. 

Since the histogram cost was one particular “histogram moving scheme”, and as we argued 
above, the “minimum permuted Li distance” is the minimum over all such schemes, we conclude 
that this Li distance is at most twice the relative earthmover distance, as desired. 

□ 


C Proof of Lemma |1] 

For convenience, we restate the lemma: 
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LemmaQ] Given two (possibly generalized) histograms g, h, a number of samples k, and a threshold 

re (0,1], 


£ (i-(!-*)*)•«(*) 

x:g{x)^0 


£ (l-(l-xt)-h(x) 

x:h(x)^ 0 


< (0.3(fc - l) + 1 )Rr(g, h) + 


Proof. We prove the inequality by considering each step of an earthmoving scheme that transforms 
g to h, and show that if in one step m probability mass is moved, at r-truncated relative earthmover 
cost r, then the sum Yl x - g ( — (1 — x) k ) ■ g(x) changes by at most (1 + 0.3 (k — 1)) • r + mkr , 
meaning that an entire earthmoving scheme to transform g into h with total cost R r (g, h) and total 
mass at most 1 changes the g term on the left hand side into the h term on the left hand side by 
changing it at most (1 + 0.3(&; — 1)) • R r (g, h) + kr. 

To prove this we first analyze the region of probability below r. By the definition of a histogram, 
m units of probability mass at probability x corresponds to a histogram entry h(x ) = and 
binomial bounds yield ^(1 — (1— x) k ) € [km(l—x^-), km], which means that when an earthmoving 
scheme moves m mass in the range x € (0, r], the expression —(1 — (1 — x ) k ) changes by at most 
km^f-r. Thus, summed over the entire earthmoving scheme, where the mass moved sums to at 
most 1, the change in Yl x - g ( x )^ofi — (1 — x ) k ) ■ g(x) from changes below probability r is at most 
k^p-r. 

To bound the remaining term, changes in Ylx-g{x)^o ^~(1 — x ) k )'9( x ) from changes in probability 
above r in the earthmoving scheme, we note that to move probability mass m from probability 
value x to y costs m\ logx — logy| in the earthmoving scheme, and changes the sum by 

We bound the ratio of these last two expressions by 1 + 0.3 (k — 1), in order to bound the total 
contribution of the portion of the earthmoving scheme above probability r by (l+0.3(fc— l))R r (g, h ), 
yielding the desired overall bound. 

We thus seek to bound the maximum change in ^ (l — (1 — relative to the change in logx 
as x changes, namely the maximum ratio of their derivatives, where we add a negative sign since 
^ (l — (1 — x) fe ) is a decreasing function. Since ^ logx = 1/x, the ratio of derivatives is 

d (1 - (1 - x) k ) _ 1 - (1 - x) fc ~ 1 ((fc ~ l)s + 1) m 

dx x x 

Consider the approximation (1 — x) fc_1 ~ e ~ x i k ~ l ), Taking logarithms of both sides, and using 
the fact that, for x < |, we have log 1 — x > — x — x 2 , we have that for x < \ the inequality 
(k — l)log(l — x) > — (k — l)(x + x 2 ); exponentiating yields (1 — x) fc_1 > e ~ x ( k ~ 1 ) . e -x 2 (k-i) > 

e -x(k- l)(f _ 

Thus for x < \ the ratio of derivatives is bounded as 

_ d (1- (l-x) fc ) 1 - (g-^lfc- 1 )(i _ x 2 (fc - l)))((fc - l)x + 1) 
dx x ~ x 

1 — e~ x ( k ~ 1 \(k — l)x + 1) e~ x ^ k ~ 1 ^ x 2 {k — l)((fe — l)x + 1) 

x x 


22 












The first term of the right hand side, after dividing by k — 1, can be reexpressed in terms of 
y = x{k — 1) as l ~ e y +1 ' , which has a global maximum less then 0.3; the second term in the right 
hand side, after the same variable substitution, equals e~ v y{y+ 1), which has a global maximum less 
than 1. Thus, for x < the absolute value of the ratio of derivatives is bounded as 0.3(fc — 1) + 1. 
For x > the right hand side of Equation |T] is \ minus some positive quantity, and is hence at 
most 2. Since 0.3(£; — 1) + 1 > 2 for any k > 5, all that remains is to checking the k = 2, 3,4 cases 
where 0.3(fc — 1) + 1 < 2 by hand to confirms that 0.3(A; — 1) + 1 is in fact a global bound. □ 

D Proof of Theorem [2] 

In this section, we prove Theorem [2 characterizing the performance of Algorithm [T] which recovers 
an accurate approximation of the histogram of the true distribution. For convenience, we restate 
Theorem [2j 


Fact [Tj There exists an absolute constant c such that for sufficiently large n and any w € [1, log ?z], 
given n independent draws from a distribution p with histogram h, with probability 1 — e _nf2{1) the 
generalized histogram hpp returned by Algorithm^ satisfies 


FI w 

n log n 


(h , Flp) < —j=- 

\/W 


The proof decomposes into three parts. In Appendix lD.il we compartmentalize the probabilistic 
portion of the proof by defining a set of conditions that are satisfied with high probability, such 
that if the samples in question satisfy the properties, then the algorithm will succeed. This section 
is analogous to the definition of a “faithful” set of samples of Definition 0 and we re-use the 
terminology of “faithful”. In Appendix ID. 21 we show that, provided the samples in question are 
“faithful”, there exists a feasible solution to the linear program defined in Algorithm [2 which 
1) has small objective function value, and 2) is very close to the true histogram from which the 
samples were drawn, in terms of T-truncated relative earthmover distance—for an appropriate 
choice of r. In Appendix ID. 31 we show that if two feasible solutions to the linear program defined in 
Algorithm [T] both have small objective function value, then they are close in tau-truncated relative 
earthmover distance. The key tool here is a Chebyshev polynomial earthmover scheme. Finally, 
in Appendix ID. 41 we put together the above pieces to prove Theorem [2j given the existence of a 
feasible point that has low-objective function value that is close to the true histogram, and the fact 
that any two solutions that both have low objective function value must be close to each other, it 
follows that the solution to the linear program that is found in Algorithm |T| must be close to the 
true histogram. 

D.l Compartmentalizing the Probabilistic Portion 

The following condition defines what it means for a set of samples drawn from a distribution to be 
“faithful” with respect to positive constants B,T> E (0,1): 

Definition 10. A set of n samples with fingerprint }~, drawn from a distribution p with histogram 
h, is said to be faithful with respect to positive constants B,V e (0,1) if the following conditions 
hold: 
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• For all i 


Ti — ^2 h(x) ■ poi{nx , i) 

x:h(x)^0 


< max ( Ff +V , n BI ' 2 + ' D '> 


For all domain elements i, letting p(i ) denote the true probability of i, the number of times i 
occurs in the samples from p differs from n ■ p{i) by at most 


max ^(n • p(i)) 2+V , n B P +v ^ . 


The “large” portion of the fingerprint F does not contain too many more samples than ex¬ 
pected: Specifically, 


Fi < n l ^ 2+v + n ^ x ■ h(x). 

i>n B +2n c x< nB f nC :h(x )>0 


The following proposition is proven via the standard “Poissonization” technique and Chernoff 
bounds. 


Proposition 2. For any constants B,V € (0,1), there is a constant a > 0 and integer uq such that 
for any n > no, a set of n samples consisting of independent draws from a distribution is “faithful” 
with respect to B,T> with probability at least 1 — e~ n . 


Proof. We first analyze the case of a Poifn )-sized sample drawn from a distribution with histogram 
h. Thus 

E[.Fj] = h(x)poi(nx,i). 

x:h(x)^0 

Additionally, the number of times each domain element occurs is independent of the number of 
times the other domain elements occur, and thus each fingerprint entry F% is the sum of independent 
random 0/1 variables, representing whether each domain element occurred exactly i times in the 
samples (i.e. contributing 1 towards F t ). By independence, Chernoff bounds apply. 

We split the analysis into two cases, according to whether E/7 7 ;] > nP. In the case that E/Fj] < 
vP , we leverage the basic Chernoff bound that if X is the sum of independent 0/1 random variables 
with E[X] < S, then for any 6 € (0,1), 


Pt[\X - E[X]\ > 5S} < 2e~ 52s/3 . 


Applied to our present setting where Fi is a sum of independent 0/1 random variables, provided 
E [Ff\ < n B , we have: 


Pr 


\Fi-E\FiW > (n B ) 2 


B\k+V 


< 2e 


(nB) 1 / 2 -® 


= 2e“ 




/3 


In the case that E/J 7 / > n , the same Chernoff bound yields 


Pr 


\Fi — E[~F/]| > E[Fi\^~ 


' E[J^ 


< 2e \nx i ] 1/2 -' D ) 3 =2e~( E[:Fi]2T ’") /3 <2e~ n2BT,/3 . 
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A union bound over the first n fingerprints shows that the probability that given a set of samples 
(consisting of Poi{n) draws), the probability that any of the fingerprint entries violate the first 

n(1) 

condition of faithful is at most n ■ 2e 3 < e n as desired. 

For the second condition of “faithful”, in analogy with the above argument, for any A < S, and 

5 € ( 0 , 1 ), 

Pr[|Po*(A) - A| > 5S] < 2e~ 52s/3 . 

Hence for x = n ■ p(i ) > n , the probability that the number of occurrences of domain element i 
differs from its expectation of n-p(i ) by at least ( n-p{i))z + ' D is bounded by 2e~^ n ' p ^ 2D l 3 < e -nn(1) . 
Similarly, in the case that x = n ■ p(i ) < n B , 

Pr[| Poi{x) — x\> n B ^ +v \ < e~ n ( \ 

For the third condition, by the Poisson tail bounds of the previous paragraph, the total aggregate 
number of occurrences of all elements with probability greater than n ^ n will differ from its 
expectation by at most n 1//2+:D , with probability 1 — e~ n ° W . Additionally, by the first condition 
of “faithful”, with probability 1 — e~ n ; no domain element i with p(i) < n will appear 
more than n B + 2 n c . Hence with probability 1 — all elements that contribute to the sum 

Yli>n B +2n c will have probability greater than n . The third condition then follows by a union 
bound over these two e -nn(1) failure probabilities. 

Thus we have shown that provided we are considering a sample size of Poi(n), the probability 
that the conditions hold is at least 1 — e -riQ<11 . To conclude, note that Pi[Poi(n) = n] > 3 - 7 ^, an d 
hence the probability that the conditions do not hold for a set of exactly n samples (namely, the 
probability that they do not hold for a set of Poi(n) samples, conditioned on the sample size being 
exactly n), is at most a factor of 3 yfn larger, and hence this probability of failure is still e~ n , as 
desired. □ 


D.2 Existence of a Good Feasible Point 

Proposition 3. Provided T is a “faithful” fingerprint derived from a distribution with histogram 
h, there exists a feasible point, (v x ), for the linear program of Algorithm]]] with objective function 
value at most 0(n^ +B+v ) such that for any r > 1 /n 3 / 2 , the r-truncated relative earthmover distance 
between the generalized histogram corresponding to (v x ) with the empirical fingerprint Fi >n B + 2 n c 

appended, and the true histogram, h, is bounded by O |max(n^ 6 ^'^,n“^ _ ^j , where the big O 
hides an absolute constant. 

Proof. Let (v x ) be defined as follows: initialize (v x ) to be identically zero. For each y < nB + nC s.t. 
h{y) > 0, increment v x by h(y)|, where x = min{x € X : x > y}. Finally, define 


m := 1 — 


E ~ Ti + 2 x • v * 

\i>nP-\-2n c x(]X 


If m > 0, increment v x by m/x for x = nB + nC . If m < 0, then arbitrarily reduce v x until a total of 
m units of mass have been removed. 
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We first argue that the r-truncated relative earthmover distance is small, and then will argue 
about the objective function value. Let h! denote the histogram obtained by appending the empiri¬ 
cal fingerprint F i>n B + 2 n c to (v x ). We construct an earthmoving scheme between h and h! as follows: 
1) for all y < nB + nC s.t. h(y) > 0, we move h(y ) • y mass to location x = min{x E X : x > y}; 2) 
for each domain element i that occurs more than nP + 2 n c times, we move p(i) mass from loca¬ 
tion p(i ) to — where X % denotes the number of occurrences of the zth domain element; 3) finally, 
whatever discrepancy remains between h and h! after the first two earthmoving phases, we move 
to probability —. Clearly this is an earthmoving scheme. For r > 1/n 3 / 2 , the r-truncated relative 

earthmover cost of the hrst phase is trivially at most log = 0(1/y/n). By the second 

condition of “faithful”, the relative earthmover cost of the second phase of the scheme is bounded 
by log(—~) = 0 (n~ B( ' 2 _:D )). To bound the cost of the third phase, note that the hrst phase 
equates the two histograms below probability vrn. By the second condition of “faithful”, after the 
second phase , there is at most 0(n~ B ^~ v ^) unmatched probability caused by the discrepancy 
between ^ and p(i) for elements observed at least nP + 2 nP times. Hence after this 2 -I b) 

discrepancy is moved to probability p-. the entirety of the remaining discrepancy lies in the prob¬ 
ability range [/—, c], where c is an upper bound on the true probability of an element that does not 
appear at least nP + 2 n c times; from the second condition of “faithful”, c < wS + 4wC ; and hence the 
total r-truncated relative earthmover distance is at most O j , as desired. 

To complete the proof of the proposition, note that by construction, (v x ) is a feasible point for 
the linear program. To see that the objective function is as claimed, note that |^poi(nx, *)| < n, 
and since we are rounding the true histogram to probabilities that are multiples of 1 /re 2 , each 
“fingerprint expectation”, X^e\'P°*( na vO ■ u x differs from Y/ x -h(x)^ oP°i( nx n) ’ h(x) by at most 
1/yfn. Together with the hrst condition of “faithful” which implies that each of the observed 
fingerprints F satisfies \Ti~ poi(rex,i)-/i(x)| < nz +T> , we conclude that the total objective 
function value is at most n B (n^ +T> + 1 /y/n) = 0(n^ +B+rD ). □ 

D.3 The Chebyshev Bump Earthmoving Scheme 

Proposition 4. Given a “ faithful” fingerprint Fi, then any pair of solutions v x ,v x to the linear 
program of Algorithm^ that both have objective function values at most 0(n^ +B+T> ) satisfy the fol¬ 
lowing: for any w € [ 1 , log re], their n ^ n -truncated relative earthmover distance R w / n \ogn[ v x,v x ] < 

0( 1/Vw). 

The proof of the above proposition relies on an explicit earthmover scheme that leverages a 
Chebyshev polynomial construction. The two key properties of the scheme are 1) the truncated 
relative earthmover cost of the scheme is small, and 2 ) given two histograms that have similar 
expected fingerprints (i.e. for all i < re , Y/ x v x poi(nx , i) « J2 x v' x poi(nx, i),) the results of applying 
the scheme to the pair of histograms will result in histograms that are very close to each other in 
truncated relative earthmover distance. We outline the construction and key propositions below. 

Definition 11. For a given re, a /3-bump earthmoving scheme is defined by a sequence of positive 
real numbers {ci}, the bump centers, and a sequence of functions {fi} : ( 0 , 1 ] —>• M such that 
E£o/i(*) = 1 f or eac h x > an d eac h function fi may be expressed as a linear combination of 
Poisson functions, fi(x) = a ijP 0 P nx ij)> such that YPjL o \ a v I < fi- 
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Given a generalized histogram h, the scheme works as follows: for each x such that h{x) / 0, 
and each integer i > 0, move xh(x ) • fi(x) units of probability mass from x to C{. We denote the 
histogram resulting from this scheme by ( c,f)(h ). 

Definition 12. A bump earthmoving scheme (c, /) is [e, T]-good if for any generalized histogram 
h the r-truncated relative earthmover distance between h and ( c,f)(h ) is at most e. 

Below we define the Chebyshev bumps to be a “third order” trigonometric construction: 

Definition 13. The Chebyshev bumps are defined in terms of n as follows. Let s = 0.2 log n. 
Define gi(y) = J2jZ- s cos(jy). Define 

92<s) = m ( 911ly ~ + 3911 2i ] + 391(9 + 2J> +91 ' (S + Ys ) ) ' 

and, for i E {1,... ,s - 1} de/ine ^(y) := y 2 (y - ^) + ff 2 (y + and $ = 92(y), and y| = 
92 ( 2 / + 7r). het tj(x) 6e f/ie linear combination of Chebyshev polynomials so that ti(cos(y)) = g\(y)- 
We thus define s + 1 functions, the “skinny bumps’’, to be Bi{x ) = ti( 1 — ^j=oP°*( xn ’ d)> f or 

i £ {0,..., s}. That is, Bi(x) is related to g\{y) by the coordinate transformation x = ^(1 —cos (y)), 
and scaling by poi(xn,j). 

Definition 14. The Chebyshev earthmoving scheme is defined in terms of n as follows: as in 
Definition WA let s = 0.2log n. For i > s + 1, define the ith bump function fi(x) = poi(nx,i — 1) 
and associated bump center Ci = . For i € {0,..., s} let fi(x) = Bi(x), and for i € {1,..., s}, 

define their associated bump centers c* = ^(1 — cos(^)), with co = c\. 

The following proposition characterizes the key properties of the Chebyshev earthmoving scheme. 
Namely, that the scheme is, in fact, an earthmoving scheme, that each bump can be expressed as a 
low-weight linear combination of Poisson functions, and that the scheme incurs a small truncated 
relative earthmover cost. 

Proposition 5. The Chebyshev earthmoving scheme of Definition \lf\ defined in terms of n, has 
the following properties: 

• For any x > 0, 

J2fi{x) = 1 , 

i>0 

hence the Chebyshev earthmoving scheme is a valid earthmoving scheme. 

• Each Bi{x) may be expressed as a ijP°i( nx ij) f or a ij satisfying 

OO 

\aij\ < 2 n 0,3 . 

j =0 


• The Chebyshev earthmoving scheme is 


0( l/y/w), 


the O notation hides an absolute constant factor. 


W 

n logn 


-good, for any w E [l,logn], where 
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The proof of the first two bullets of the proposition closely follow the arguments in [29]. For 
the final bullet point, the intuition of the proof is the following: the ?'th bump Bi, with center 
Ci = ^ (1 — cos(*7r/s)) ~ i 2 ^- s has a width of O(^), and Bi(x ) decays rapidly (as the fourth 
power) away from its center, Cj. Specifically, F>j(cj ± ^|) < 0(l/a 4 ). Hence, at worst, the cost of 
the earthmoving scheme will be dominated by the cost of moving the mass around the smallest c* 
that exceeds the truncation parameter w/n log n. Such a bump will have width O(^) = 0( n ^ n ), 
which will incur a per-unit mass relative earthmover cost of 0(y/l/w). 

For completeness, we give a complete proof of Proposition [5l with the three parts split into 
distinct lemmas: 


Lemma 8. For any x 

s . oo 

22 92(x-\ - : ) = 1 , and 22 fi(x) = 1 . 

i=— s+1 2=0 

Proof. 52 ( 2 /) is a linear combination of cosines at integer frequencies j, for j = 0 , ,s, shifted by 
±7 t/2s and ±37r/s2. Since J2i=- s +i 9z( x + ir) sums these cosines over all possible multiples of tt/s, 
we note that all but the frequency 0 terms will cancel. The cos(Oy) = 1 term will show up once in 
each gi term, and thus l + 3 + 3-|-l = 8 times in each 52 term, and thus 8 • 2s times in the sum in 
question. Together with the normalizing factor of 16s, the total sum is thus 1, as claimed. 

For the second part of the claim, 


E /<( 


i x = 


xn 


i=o 


2s 


-1 + 


7TJ, 


s— 1 


22 92 (cos 1 

\J=-s+1 

S— 1 

1 • ^^poz(xn, j) + ^^poi(xn, j) = 1. 

j =0 j>s 


"22poi(xn,j ) + 22P oi ( xn ’i) 
j =0 j>s 


□ 

We now show that each Chebyshev bump may be expressed as a low-weight linear combination 
of Poisson functions. 

Lemma 9. Each Bi(x ) may be expressed as a ijP°i(nx,j) for aij satisfying 

OO 

22 \ a ij I ^ 2n 0 ' 3 . 

3 =0 

Proof. Consider decomposing g\(y) into a linear combination of cos (iy), for t € {0,... ,s}. Since 
cos (—iy) = cos (£y), gi(y) consists of one copy of cos(sy), two copies of cos (£y) for each i between 
0 and s, and one copy of cos(Oy); 52 ( 2 /) consists of (-^= times) 8 copies of different ^i(y)’s, with 
some shifted so as to introduce sine components, but these sine components are canceled out in the 
formation of 223 ( 2 /), which is a symmetric function for each i. Thus since each <73 contains at most 
two < 72 ’s, each 2 / 3 ( 2 /) ma Y be regarded as a linear combination ^| =0 cos(ly)bu with the coefficients 
bounded as \ba\ < 
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Since t,; was defined so that fj(cos(y)) = g\(y) = X)l=o cos (£y)ba, by the definition of Chebyshev 
polynomials we have U(z) = Ye=oTi{ z )bu- Thus the bumps are expressed as 


*(*)= £r«a-£ 


a =o 



3-1 


5 ~2poi(xn,j ) 
j=o 


We further express each Chebyshev polynomial via its coefficients as Tj?( 1 — (1 — 

|j) m and then expand each term via binomial expansion as (1 — |f) m = Ylq= o( — lf) 9 Cq) to yield 


s i m s —1 

= Z Z Z Z (-^7 

£=0 m =0 q =0 j =0 


xn\<? (m 

~27) 


bit poi(xn,j). 


J, ^ —v— 

finally lets us express Bj as a linear combination of Poisson functions, for all i G {0,..., s}: 


We note that in general we can reexpress x q poi(xn, j) = x q xJ ' nJ q - = poi(xn,j + q) , which 


s £ m s —1 

£*(*) = ZZZZ' 

£=0 m =0 5=0 j =0 


2 s 


m\ (j + q )! 


Q 


b ie poi(xn,j + q). 


It remains to bound the sum of the absolute values of the coefficients of the Poisson functions. 
That is, by the triangle inequality, it is sufficient to show that 


s t m s —1 

££££ 

£=0 m =0 5=0 j=0 


1 

'Ys 


m\ (j + q )! 


bu 


< 2 n 


0.3 


We take the sum over j first: the general fact that Ym =o (Z*) = Ot+i ) iurplies that 
YjZl-jT- = Ej=o (Z) g! = ?!(Si) = 5+1 ( 3 -?jp and further, since q < m < £ < s we have 
s + q < 2 s which implies that this final expression is bounded as jf~[yT = s 7 Pi Zw — s ' ( 2s ) 9 - 


Thus we have 


s £ m s —1 

££££ 

£=0 m =0 5=0 j =0 


1 \ q (m\ (j + q)l 


2s 


j! 


-bi£ 


< 


s i m 

£££ 

£=0 m =0 q —0 




m 


bu 


1 Z im Z 


£=0 


m— 0 


Chebyshev polynomials have coefficients whose signs repeat in the pattern (+,0, — ,0), and 
thus we can evaluate the innermost sum exactly as |T7(2r) |, for r = \f— T. Since we bounded 
\bie\ < | above, the quantity to be bounded is now s Y,£=o f |T^(2r)|. Since the explicit expression for 
Chebyshev polynomials yields |TJ»(2r)| = \ [(2 — \fh) £ + (2 + V5) £ ] and since |2— \fh\ £ = (2+\/5)~^ 

we finally bound s^=o § \ T i( 2r )\ < 1 + EL-s ( 2 + < 1 + 2 +VE^-i ’ ( 2 + ^) S < 2 ‘ ( 2 + v / 5) s < 

2 • A; 0 - 3 , as desired, since s = 0.2 log n and log(2 + \/5) < 1.5 and 0.2 • 1.5 = 0.3. □ 

The following lemma quantifies the “skinnyness” of the Chebyshev bumps, which is the main 
component in the proof of the quality of the scheme (the third bullet in Proposition [5]). 
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7 

Lemma 10. \g 2 (y)\ < ^ for y E [—7r, 7r] \ (-37t/s, 37t/s), and \g 2 (y)\ < 1/2 everywhere. 

Proof. Since gi(y) = cos J5 = sin(sy) cot(y/2), and since sin(a + 7r) = — sin(a), we have the 

following: 


52 (y) = 


1 ( . 37T. . 7T . 7T . , 37T. 

Ifo ( 9lfe - 27* + 3 » l( » - + 91 < 9 + 2j) 

ik ( sinfes+T/2) ( cot( l -1/ - 3cot( l - 


Note that (cot(| — ||) — 3cot(| — |j) + 3cot(| + — cot (If + ||)) is a discrete approximation 
to (7t/2s) 3 times the third derivative of the cotangent function evaluated at y/2. Thus it is bounded 
in magnitude by (7t/2s) 3 times the maximum magnitude of 4^ cot(x) in the range x E [| — | + 

||]. Since the magnitude of this third derivative is decreasing for x E (0, it), we can simply evaluate 
the magnitude of this derivative at | —We thus have cot (x) = , whose magnitude 

is at most ^ or x ^ (0, 7r] . For y E [37t/s,7t], we trivially have that | — and thus we 

have the following bound: 


,y 37r. „ .7/ 

cot(-) — 3 cot( — 

v 2 4s v 2 


7r , „ ,y 7r , 

—) + 3 cot(—|- 

4s' v 2 4s' 


,y 37r.. 

cot( 2 + u 11 - 


/ 7T \ 3 6 127T 7 

V2s/ (y/ 2tt) 4 ~ y 4 s 3 


Since 52 ( 5 ) is a symmetric function, the same bound holds for y E [— ir, — 2>tt / s\. Thus ( 52 ( 5)1 < 
16s 2 y 4 s a < for d ^ [ —7r i 7r ] \ (~ 37r/s, 3n/s). To conclude, note that 52(5) attains a global 
maximum at y = 0, with 52 ( 0 ) = (6cot(7r/4s) — 2cot(37r/4s)) < < 1/2. □ 

We now prove the final bullet point of Proposition [5] 


Lemma 11. The Chebyshev earthmoving scheme is 


0(l/y/w), 


where the O notation hides an absolute constant factor. 


n log 71 


-good, for any w E [1, log n], 


Proof. We split this proof into two parts: first we will consider the cost of the portion of the scheme 
associated with all but the first s + 1 bumps, and then we consider the cost of the skinny bumps /,; 
with i E {0,..., s}. 

For the first part, we consider the cost of bumps fi for i > s +1; that is the relative earthmover 
cost of moving poi(xn , i) mass from x to summed over i> s. By definition of relative earthmover 
distance, the cost of moving mass from x to p- is | log ^|, which, since logy < y — 1, we bound by 
/p — 1 when i < xn and — 1 otherwise. We thus split the sum into two parts. 

For i > [xn] we have poi(xn, i )(—— 1) = poifxn, i — 1) — poifxn, i). This expression telescopes 
when summed over i > maxjs, [xn]} to yield poi(xn, max{s, [xn]} — 1) = 

For i < [xn] — 1 we have, since i > s, that poi(xn,i)(™ — 1) < poi(xn,i)(( 1 + l)4pt- — 1) = 
(1 + -)poi(xn,i + 1) — poi(xn,i). The 1 term sums to at most p, and the rest telescopes to 
poifxn , [xn]) — poi(xn, s) = O(-^j). Thus in total, fi for i > s + 1 contributes O(-^j) to the relative 
earthmover cost, per unit of weight moved. 
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We now turn to the skinny bumps fi(x) for i < s. The simplest case is when x is outside 
the region that corresponds to the cosine of a real number — that is, when xn > 4s. It is 
straightforward to show that fi{x) is very small in this region. We note the general expression 
for Chebyshev polynomials: Tj{x) = ^ (x — Vx 2 — l)- 7 + (x + Vx 2 — l)- 7 , whose magnitude we 
bound by |2xp'. Further, since 2x < ^e x , we bound this by (f pe^- 7 , which we apply when \x\ > 1. 
Recall the definition fi{x) = £j(l — ff) X^j=oP 0? '( xn )i)) where t % is the polynomial dehned so that 
ti(cos(y)) = g l 3 (y), that is, tj. is a linear combination of Chebyshev polynomials of degree at most s 
and with coefficients summing in magnitude to at most 2, as was shown in the proof of Lemma 0 
Since xn > s, we may bound Ylj=oP°4 xn , j) < s-poi(xn, s). Further, since z < e z ~ l for all z, letting 


z = yields x < 4s • 1 

4 s s s < As p —3xn/4 
e s. e 3xn/4 c! — t: c- 


, from which we may bound poi(xn, s) = 


( xn) s e 


< 


. , xn -i . „ 

-(4s • e 4s ) = 


We combine this with the above bound on the magnitude of Chebyshev 
polynomials, Tj(z) < (|pel z h < (|) s el z l s , where z = (1 — yields Tj(z) < (p ) s e~. Thus 
fi(x ) < poly(s)4 s e~ 3xn ^ 4: (p) s e^~ = poly(s)(p) s e^~. Since ™ > s in this case, fi is exponentially 
small in both x and s; the total cost of this earthmoving scheme, per unit of mass above — 
is obtained by multiplying this by the logarithmic relative distance the mass has to move, and 
summing over the s + 1 values of i < s, and thus remains exponentially small, and is thus trivially 
bounded by 0(-f=). 

To bound the cost in the remaining case, when xn < 4s and i < s, we work with the trigono¬ 
metric functions g l 3 , instead of L directly. Since mass may be moved freely below probability n ^ n , 
we may assume that all the mass below this value is located at probability exactly n ^ gn - 

For y € (0,vr], we seek to bound the per-unit-mass relative earthmover cost of, for each i > 0, 
moving g 3 (y) mass from — (1 — cos (y)) to q. By the above comments, it suffices to consider 
y E [0(- y ^),7r]. This contribution is at most 


1 03(y) f 1 °g( 1 - cos (y)) - iogl 1 - cos(—)) 

t=0 v s 

We analyze this expression by first showing that for any x,x' € (0,7r], 

|log(1 — cos(x)) — log(l — cos(x / ))| < 2| logx — log x' I. 

Indeed, this holds because the derivative of log(l — cos{x)) is positive, and strictly less than 
the derivative of 21ogx; this can be seen by noting that the respective derivatives are and 

|, and we claim that the second expression is always greater. To compare the two expressions, 
cross-multiply and take the difference, to yield ysiny — 2 + 2cos y, which we show is always at 
most 0 by noting that it is 0 when y = 0 and has derivative y cos y — sin y, which is negative since 
y < tan y. Thus we have that | log(l — cos (y)) — log(l — cos(^))| < 2\ log y — log ^|; we use this 
bound in all but the last step of the analysis. Additionally, we ignore the ^2j=lpoi(xn, j) term as 
it is always at most 1. 

We will now show that 

\g 3 (y)(logy- log - )| + Y \9\(y)4°gy ~ log — )| = 0(—), 
s 4—' s sy 



where the first term is the contribution from /o,co- For i such that y E 3 ^ , ( j+3 ) 7r ) ) by the 
second bounds on l^l in the statement of LemmaCOIl g\{y) < 1, and for each of the at most 6 such 
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i, |(logy — log max l 1 ’ ; l 7r )| < _L ; to yield a contribution of O(^). For the contribution from i such 

that y < or y > ; the first bound of Lemma HOl yields |^3 (y) | = Roughly, 

the bound will follow from noting that this sum of inverse fourth powers is dominated by the first 
few terms. Formally, we split up our sum over i £ [s] \ — 3, ^ + 3] into two parts according to 

whether i > ys/ir. 


E l - n ITT. . ^—a 7r / . 

7 - —4 logy -log— < > Ws - -4 log* -log — 

V. , „ (y s - m ) S (*- - i) 4 7T 


< 7T 


/ 

•7 ID 


1 n , V s ■ 

775 - ^4 (log u; — log 

=a£+2 tt 


(2) 


Since the antiderivative of (log w — log a) with respect to w is 


—2 w(w 2 — 3 wa + 3a 2 ) log w + 2(w — a) 3 log(tc — a) + a(2ui 2 — 5 wa + 3a 2 + 2a 2 log a) 

6(u; — a) 3 a 3 

the quantity in Equation [2] is equal to the above expression evaluated with a = and w = a + 2, 
to yield 

0(1) - log V - + log (2 + V -) = 0(1). 
ys it tt ys 

A nearly identical argument applies to the portion of the sum for i < ^ + 3, yielding the same 
asymptotic bound of O(^). As it suffices to consider y > O(^), this bounds the total per-unit 
mass -r 1 -—truncated relative earthmover cost as 0( -4=), as desired. 

n log n \ y/w ' ’ 

□ 


D.4 Proof of Theorem \ 2 \ 


We now assemble the key propositions from the above sections to complete our proof of Theorem [2j 

Proposition [2] guarantees that with high probability, the samples will be “faithful”. For the 
remainder of the proof, we will assume that we are working with a faithful set of n independent 
draws from a distribution with true histogram h. Proposition [2] guarantees that there exists a 
feasible point ( v x ) for the linear program of Algorithm [T] with objective function at most 0(n^ +l3+c ), 
such that if the empirical fingerprint above probability n +2n is appended, the resulting histogram 
hi satisfies R r (h,hi) < O |max(n' 6 ^"®^,n~( K ” c ))j , for any r > 1 /n 3 / 2 . 

Let h 2 denote the histogram resulting from Algorithm [l] Hence the portion of h 2 below proba¬ 
bility nB +^ nC corresponds to a feasible point of the linear program with objective function bounded 
by 0 (n 2 +B+c ). Additionally, h\{x) and h 2 {x) are identical for all x > nB + nC , as, by construction, 
they are both zero for all x € ( n , n + ^ n ], and are both equal to the empirical distribution 
of the samples above this region. We will now leverage the Chebyshev earthmoving scheme, via 
Proposition [5] to argue that for any w € [l,lognl, R m (hi. ho) < O (-)=), and hence by the 

triangle inequality, R. 


^_(h,h 2 )<0(^). 


log r 


To leverage the Chebyshev earthmoving scheme, recall that the earthmoving scheme that moves 
all the probability mass of a histogram to a discrete set of “bump centers” (q), such that the earth 
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moving scheme incurs a small truncated relative earthmover distance, and also has the property 
that when applied to any histogram g, the amount of probability mass that ends up at each bump 
center, q is given as Yj>o Ei-jW/o P°K nx ,j)xg(x), for some set of coefficients a t .j satisfying 
for all i. Y 


j> o 


a 


< 2n 0 ' 3 . 


Consider the results of applying the Chebyshev earthmoving scheme to histograms hi and h 2 . 
We first argue that the discrepancy in the amount of probability mass that results at the ith bump 
center will be negligible for any i > nP + 2rf. Indeed, since h\ and h 2 are identical above probability 
n and Yi>n B +- 2 C P°i(Aj *) = e ~ n for A < rP + n c , the discrepancy in the mass at all bump 
centers q for i >> nP + 2 n c is trivially bounded by o(l/n). 

We now address the discrepancy in the mass at the bump centers c* for i < rP + 2 rP. For any 
such i the discrepancy is bounded by the following quantity: 


poi (nx,j)x (h\(x) — h 2 (x)) 

j> 0 x:h(x)^ 0 


E E -poi (nx,j + 1) (hi(x) - h 2 (x)) 

j >0 x:h(x )^0 


< 


< 


< 


' n 

j> 1 


Y poi(nx,j) (hi(x) - h 2 (x)) 

x:h(x)^0 
n e + An c 

J 


)(!/n)+ Y a hj -1 


n 


3 = 1 

) 3 (rP + 4n c ) 2 


n 


Y P°K nx J) ( h i( x ) - h 2 (x)) 

x\h(x)^ 0 


n 


0(n 


i+B+Cx 


0(n' 


0.3+4+3B+C—1 \ 


Where, in the third line, we leveraged the bound Yj \ a i,j\ — n °’ 3 an d the bound of 0(n^ +l3+c ) on 
the linear program objective function corresponding to h\ and h 2 , which measures the discrepancies 
between Y x poi (nx,j)h.(x) and the corresponding fingerprint entries. Note that the entirety of this 
discrepancy can be trivially equalized at a relative earthmover cost of 


0(n°- 3- 4 +3B+C ~ 1 log(n)), 


by, for example, moving this discrepancy to probability value 1. To complete the proof, by the 
triangle inequality we have that for any w € [l,logn], letting g\ and g 2 denote the respective 
results of applying the Chebyshev earthmoving scheme to histograms h\ and h 2 . we have the 
following: 


n log n 


(h,h 2 ) 


< R™_(h,hi) + R™_(hi,gi) + R*>_(gi, g 2 ) + R_m_(g 2 , h 2 ) 

n log n n log n n log n n log n 

< O |ma x(n~ s ^~ T> ') } n~P~ c ^)J + 0(1/y/w) + O(n 0 ' 3+ 2 +3B+c_1 log(n)) + 0(1/y/w) 

< 0(l/yfw). 


E Rounding a Generalized Histogram 

Algorithm |Tj returns a generalized histogram. Recall that generalized histograms are histograms 
but without the condition that their values are integers, and thus may not correspond to actual 
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distributions—whose histogram entries are always integral. While a generalized distribution suffices 
to establish Theorem [lj we observe that it is possible to round a generalized histogram without 
significantly altering it, in truncated relative earthmover distance. The following algorithm and 
lemma characterizing its performance show one way to round the generalized histogram to obtain a 
histogram that is close in truncated relative earthmover distance. This, together with Theorem (2] 
establishes Corollary [U 


Algorithm 3. Round to Histogram 
Input: Generalized histogram g. 

Output: Histogram h. 

• Initialize h to consist of the integral elements of g. 

• For each integer j > 0: 

— Let Xji , Xj 2 , ■ ■ ■, Xjt be the elements of the support of g that lie in the range [ 2 ^+!). 2“- J ] 
and that have non-integral histogram entries; let m := X^i=i x ji9( x ji) be the total mass 
represented; initialize histogram h' to be identically 0 and set variable diff := 0. 

— For i = 

* If diff < 0 set h'(xji ) = \g( x ji)], otherwise, if diff > 0 set h'(xji) = \_g( x ji)\- 

* Increment diff by Xji ( h'(xji ) — g{xji)). 

— For each i e 1 increment x ji ) by h!{xji). 


Lemma 12. Let h be the output of running Algorithm 0 on generalized histogram g. The following 
conditions hold: 

• For all x, h(x) € NU {0}, and J2 x -h(x)f=o x h( x ) = hence h is a histogram of a distribution. 

• Ro(h,g) < 20a where a := max(x : g(x) 0 NU {0}). 

Proof. For each stage j of Algorithm 03 the algorithm goes through each of the histogram entries 
g(xji) rounding them up or down to corresponding values h'(xji) and storing the cumulative dif¬ 
ference in probability mass in the variable diff. Thus if this region of g initially had probability 
mass m, then h! will have probability mass m + diff. We bound this by noting that since the first 
element of each stage is always rounded up, and 2 _( b+ 1 ) i s the smallest possible coordinate in this 
stage, the mass of h 1 , namely m + diff, is thus always at least 2 _ b/+ 1 ). Since each element of h! is 
scaled by before being added to h, the total mass contributed by stage j to h is exactly m, 

meaning that each stage of rounding is “mass-preserving”. 

Denoting by gj the portion of g considered in stage j, and denoting by hj this stage’s contribution 
to h, we now seek to bound R(hj,gj). 

Recall the cumulative distribution, which for any distribution over the reals, and any number y, 
is the total amount of probability mass in the distribution between 0 and y. Given a generalized his¬ 
togram g, we can define its (generalized) cumulative distribution by c(g)(x) := Yl x <yg(x)^o x d( x )- 
We note that at each stage j of Algorithm [3] and in each iteration i of the inner loop, the variable 
diff equals the difference between the cumulative distributions of h' and gj at xji , and hence 
also on the region immediately to the right of Xji. Further, we note that at iteration i, \diff\ 
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is bounded by Xji since at each iteration, if diff is positive it will decrease and if it is nega¬ 
tive it will increase, and since h'fxjf) is a rounded version of g(xji), diff will be changed by 
Xji(h!(xji) — g{xji )) which has magnitude at most Xji . Combining these two observations yields 
that for all x, \c(h')(x) — c(gj)(x)\ < x. 

To bound the relative earthmover distance we note that for distributions over the reals, the 
earthmover distance between two distributions can be expressed as the integral of the absolute 
value of the difference between their cumulative distributions; since relative earthmover distance 
can be related to the standard earthmover distance by changing each x value to log x, the change 
of variables theorem gives us that R(a,b) = f ^\c(b)(x) — c(a)(x)\ dx. We can thus use the bound 
from the previous paragraph in this equation after one modification: since h! has total probability 
mass m + diff, its relative earthmover distance to gj with probability mass m is undefined, and we 
thus define h" to be h! with the modification that we subtract diff probability mass from location 
2~ J (it does not matter to this formalism if diff is negative, or if this makes h''(2~ J ) negative). 
We thus have that R(h" ,gj) = f^-u+ 1 ) y| c{h!)(x) — c(gj)(x) \ dx < / 2 2 _ ( j+1 ) \xdx = 

We now bound the relative earthmover distance from h" to hj via the following two-part earth- 
moving scheme: all of the mass in h" that comes from h! (specifically, all the mass except the 
—diff mass added at 2~ 3 ) is moved to a fraction of its original location, at a relative 

earthmover cost (m + diff) ■ | log |: the remaining —diff mass is moved wherever needed, 

involving changing its location by a factor as much as 2 • max{} at a relative earth- 
mover cost of at most \diff\ ■ (log2 + | log [)• Thus our total bound on R(gj,hj), by the 

triangle inequality, is 2~^ +v > + (m + diff) ■ \ log , m+ ^ :// 1 + \diff\ ■ (log 2 + | log which we 

use when m > 2 ~ J , in conjunction with the two bounds derived above, that \diff\ < 2 _J and that 
m + diff > yielding a total bound on the earthmover distance of 5 • 2~i for the jth stage 

when m > 2~ ] . When m < 2~ J we note directly that m mass is being moved a relative distance of 
at most 2 • max{ , m+ ^^ } at a cost of m ■ (log 2 + | log I) which we again bound by 

5 • 2~ 3 . Thus, summing over all j > [| log 2 a|J, yields a bound of 20a. □ 


35 



